[PYTHON MATRIX-SIG] numexprmodule.c

Hinsen Konrad hinsenk@ere.umontreal.ca
Fri, 19 Jan 1996 16:25:48 -0500


/* Numexpr objects */

#include "allobjects.h"
#include <math.h>

/* Type of stack entries */

typedef union { long i; double f; complex c; } stack_entry;

/* Error object */

static PyObject *ErrorObject;

/* Object definition */

typedef struct {
	PyObject_HEAD
	int   nargs;		/* number of parameters */
	char *argtypes;		/* types of parameters */
	stack_entry *args;	/* array to store actual arguments */
	int   nconstants;	/* number of constants */
	stack_entry *constants;	/* constant array */
	int   codesize;		/* number of code words */
	int  *code;		/* compiled expression */
	int   stacksize;	/* stack size needed for evaluation */
	stack_entry *stack;	/* execution stack */
	char  type;		/* return type of the expression */
} PyNumExprObject;

staticforward PyTypeObject PyNumExpr_Type;

#define PyNumExpr_Check(v)  	((v)->ob_type == &PyNumExpr_Type)

/* Type codes */

#define INTEGER		1
#define FLOAT		2
#define COMPLEX		3
#define NTYPES		4

/* Opcodes used in compiled code */

#define GET_ARG		1
#define GET_CONSTANT	2
#define ADD		10
#define ADD_INT		11
#define ADD_FLOAT	12
#define ADD_COMPLEX	13
#define SUB		15
#define SUB_INT		16
#define SUB_FLOAT	17
#define SUB_COMPLEX	18
#define MUL		20
#define MUL_INT		21
#define MUL_FLOAT	22
#define MUL_COMPLEX	23
#define DIV		25
#define DIV_INT		26
#define DIV_FLOAT	27
#define DIV_COMPLEX	28
#define POW		30
#define NEG		40
#define NEG_INT		41
#define NEG_FLOAT	42
#define NEG_COMPLEX	43
#define ABS		45
#define ABS_INT		46
#define ABS_FLOAT	47
#define ABS_COMPLEX	48
#define TYPE_CAST	50
#define	INTEGER_FLOAT	59
#define	INTEGER_COMPLEX	63
#define FLOAT_INTEGER	56
#define FLOAT_COMPLEX	64
#define COMPLEX_INTEGER 57
#define COMPLEX_FLOAT   61

#define ACOS_FLOAT	100
#define ACOS_COMPLEX	101
#define ACOSH_FLOAT	102
#define ACOSH_COMPLEX	103
#define ASIN_FLOAT	104
#define ASIN_COMPLEX	105
#define ASINH_FLOAT	106
#define ASINH_COMPLEX	107
#define ATAN_FLOAT	108
#define ATAN_COMPLEX	109
#define ATANH_FLOAT	110
#define ATANH_COMPLEX	111
#define COS_FLOAT	112
#define COS_COMPLEX	113
#define COSH_FLOAT	114
#define COSH_COMPLEX	115
#define EXP_FLOAT	116
#define EXP_COMPLEX	117
#define LOG_FLOAT	118
#define LOG_COMPLEX	119
#define LOG10_FLOAT	120
#define LOG10_COMPLEX	121
#define SIN_FLOAT	122
#define SIN_COMPLEX	123
#define SINH_FLOAT	124
#define SINH_COMPLEX	125
#define SQRT_FLOAT	126
#define SQRT_COMPLEX	127
#define TAN_FLOAT	128
#define TAN_COMPLEX	129
#define TANH_FLOAT	130
#define TANH_COMPLEX	131


/* Utility functions */

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


/* Allocation and deallocation of numexpr objects */

static PyNumExprObject *
new_numexprobject()
{
	PyNumExprObject *self;
	self = NEWOBJ(PyNumExprObject, &PyNumExpr_Type);
	if (self == NULL)
		return NULL;
	self->nargs = 0;
	self->argtypes = NULL;
	self->args = NULL;
	self->nconstants = 0;
	self->constants = NULL;
	self->codesize = 0;
	self->code = NULL;
	self->stacksize = 0;
	self->stack = NULL;
	self->type = 0;
	return self;
}

static void
numexpr_dealloc(self)
	PyNumExprObject *self;
{
	PyMem_XDEL(self->argtypes);
	PyMem_XDEL(self->args);
	PyMem_XDEL(self->constants);
	PyMem_XDEL(self->code);
	PyMem_XDEL(self->stack);
	PyMem_DEL(self);
}

/* Attribute access. The only attributes are methods. */

struct PyMethodDef numexpr_methods[];

static PyObject *
numexpr_getattr(self, name)
	PyNumExprObject *self;
	char *name;
{
	return Py_FindMethod(numexpr_methods, (PyObject *)self, name);
}

/* Printed representation */

static PyObject *
numexpr_repr(self)
	PyNumExprObject *self;
{
	char buf[100];
	sprintf(buf, "<compiled expression (%s)>",
			self->argtypes);
	return PyString_FromString(buf);
}

/* Evaluation of a compiled expression */

extern complex c_sqrt();
extern complex c_acos();
extern complex c_acosh();
extern complex c_asin();
extern complex c_asinh();
extern complex c_atan();
extern complex c_atanh();
extern complex c_cos();
extern complex c_cosh();
extern complex c_exp();
extern complex c_log();
extern complex c_log10();
extern complex c_sin();
extern complex c_sinh();
extern complex c_sqrt();
extern complex c_tan();
extern complex c_tanh();

static PyObject *
numexpr_call(self, args)
	PyNumExprObject *self;
	PyObject *args;
{
	stack_entry *stack = self->stack;
	int *cp, *end;
	char text[50];
	stack_entry *sp;
	PyObject *arg;
	char argtype[2];
	PyObject *result;
	int i;

	if (PyObject_Length(args) != self->nargs) {
		PyErr_SetString(TypeError,
			"wrong number of arguments in compiled expression");
		return NULL;
	}
	argtype[1] = '\0';
	for (i = 0; i < self->nargs; i++)
	{
		arg = PySequence_GetItem(args, i);
		argtype[0] = self->argtypes[i];
		if (argtype[0] == ' ')
			PyArg_Parse(arg, "O", &arg);
		else
			PyArg_Parse(arg, argtype, self->args+i);
	}

	end = self->code + self->codesize;
	cp = self->code;
	sp = stack-1;
	while (cp < end) {
		switch (*cp) {

		case GET_ARG:
			cp++;
			*++sp = self->args[*cp];
			break;
		case GET_CONSTANT:
			cp++;
			*++sp = self->constants[*cp];
			break;

		case ADD_INT:
			sp[-1].i += sp[0].i;
			sp--;
			break;
		case ADD_FLOAT:
			sp[-1].f += sp[0].f;
			sp--;
			break;
		case ADD_COMPLEX:
			sp[-1].c.real += sp[0].c.real;
			sp[-1].c.imag += sp[0].c.imag;
			sp--;
			break;

		case SUB_INT:
			sp[-1].i -= sp[0].i;
			sp--;
			break;
		case SUB_FLOAT:
			sp[-1].f -= sp[0].f;
			sp--;
			break;
		case SUB_COMPLEX:
			sp[-1].c.real -= sp[0].c.real;
			sp[-1].c.imag -= sp[0].c.imag;
			sp--;
			break;

		case MUL_INT:
			sp[-1].i *= sp[0].i;
			sp--;
			break;
		case MUL_FLOAT:
			sp[-1].f *= sp[0].f;
			sp--;
			break;
		case MUL_COMPLEX:
			sp[-1].c = c_prod(sp[-1].c, sp[0].c);
			sp--;
			break;

		case DIV_INT:
			sp[-1].i /= sp[0].i;
			sp--;
			break;
		case DIV_FLOAT:
			sp[-1].f /= sp[0].f;
			sp--;
			break;
		case DIV_COMPLEX:
			sp[-1].c = c_quot(sp[-1].c, sp[0].c);
			sp--;
			break;

		case INTEGER_FLOAT:
			sp->f = sp->i;
			break;
		case INTEGER_COMPLEX:
			sp->c.real = sp->i;
			sp->c.imag = 0.;
			break;
		case FLOAT_INTEGER:
			sp->i = sp->f;
			break;
		case FLOAT_COMPLEX:
			sp->c.real = sp->f;
			sp->c.imag = 0.;
			break;
		case COMPLEX_INTEGER:
			sp->i = sp->c.real;
			break;
		case COMPLEX_FLOAT:
			sp->f = sp->c.real;
			break;

		case NEG_INT:
			sp->i = -sp->i;
			break;
		case NEG_FLOAT:
			sp->f = -sp->f;
			break;
		case NEG_COMPLEX:
			sp->c.real = -sp->c.real;
			sp->c.imag = -sp->c.imag;
			break;

		case ABS_INT:
			if (sp->i < 0)
				sp->i = -sp->i;
			break;
		case ABS_FLOAT:
			if (sp->f < 0)
				sp->f = -sp->f;
			break;
		case ABS_COMPLEX:
			sp->f = hypot(sp->c.real, sp->c.imag);
			break;

		case ACOS_FLOAT:
			sp->f = acos(sp->f);
			break;
		case ACOS_COMPLEX:
			sp->c = c_acos(sp->c);
			break;
		case ACOSH_FLOAT:
			sp->f = acosh(sp->f);
			break;
		case ACOSH_COMPLEX:
			sp->c = c_acosh(sp->c);
			break;
		case ASIN_FLOAT:
			sp->f = asin(sp->f);
			break;
		case ASIN_COMPLEX:
			sp->c = c_asin(sp->c);
			break;
		case ASINH_FLOAT:
			sp->f = asinh(sp->f);
			break;
		case ASINH_COMPLEX:
			sp->c = c_asinh(sp->c);
			break;
		case ATAN_FLOAT:
			sp->f = atan(sp->f);
			break;
		case ATAN_COMPLEX:
			sp->c = c_atan(sp->c);
			break;
		case ATANH_FLOAT:
			sp->f = atanh(sp->f);
			break;
		case ATANH_COMPLEX:
			sp->c = c_atanh(sp->c);
			break;
		case COS_FLOAT:
			sp->f = cos(sp->f);
			break;
		case COS_COMPLEX:
			sp->c = c_cos(sp->c);
			break;
		case COSH_FLOAT:
			sp->f = cosh(sp->f);
			break;
		case COSH_COMPLEX:
			sp->c = c_cosh(sp->c);
			break;
		case EXP_FLOAT:
			sp->f = exp(sp->f);
			break;
		case EXP_COMPLEX:
			sp->c = c_exp(sp->c);
			break;
		case LOG_FLOAT:
			sp->f = log(sp->f);
			break;
		case LOG_COMPLEX:
			sp->c = c_log(sp->c);
			break;
		case LOG10_FLOAT:
			sp->f = log10(sp->f);
			break;
		case LOG10_COMPLEX:
			sp->c = c_log10(sp->c);
			break;
		case SIN_FLOAT:
			sp->f = sin(sp->f);
			break;
		case SIN_COMPLEX:
			sp->c = c_sin(sp->c);
			break;
		case SINH_FLOAT:
			sp->f = sinh(sp->f);
			break;
		case SINH_COMPLEX:
			sp->c = c_sinh(sp->c);
			break;
		case SQRT_FLOAT:
			sp->f = sqrt(sp->f);
			break;
		case SQRT_COMPLEX:
			sp->c = c_sqrt(sp->c);
			break;
		case TAN_FLOAT:
			sp->f = tan(sp->f);
			break;
		case TAN_COMPLEX:
			sp->c = c_tan(sp->c);
			break;
		case TANH_FLOAT:
			sp->f = tanh(sp->f);
			break;
		case TANH_COMPLEX:
			sp->c = c_tanh(sp->c);
			break;

		default:
			sprintf(text, "unknown opcode %d", *cp);
			PyErr_SetString(ErrorObject, text);
			return NULL;
			break;
		}
		cp++;
       	}

	if (sp != stack) {
		PyErr_SetString(ErrorObject, "stack error");
		return NULL;
	}
	switch(self->type) {
	case INTEGER:
		result = (PyObject *)PyInt_FromLong(sp->i);
		break;
	case FLOAT:
		result = (PyObject *)PyFloat_FromDouble(sp->f);
		break;
	case COMPLEX:
		result = (PyObject *)PyComplex_FromCComplex(sp->c);
		break;
	default:
		Py_INCREF(None);
		result = None;
	}
	return result;
}

/* Arithmetic */

static char
numexpr_commontype(type1, type2)
	char type1;
	char type2;
{
	return max(type1, type2);
}

static PyNumExprObject *
numexpr_binary_op(a, b, type, op)
	PyNumExprObject *a;
	PyNumExprObject *b;
	char type;
	int op;
{
	PyNumExprObject *r = new_numexprobject();
	int error = 0;
	int i, j;
	int *cp;

	if (r != NULL) {
		r->type = type;
		r->nargs = max(a->nargs, b->nargs);
		r->nconstants = a->nconstants + b->nconstants;
		r->stacksize = max(a->stacksize, b->stacksize) + 1;
		r->codesize = a->codesize + b->codesize
			+ (a->type != type) + (b->type != type) + 1;
		if ((r->argtypes = PyMem_NEW(char, r->nargs+1)) != NULL) {
			for (i = 0; i < min(a->nargs, b->nargs); i++) {
				if (a->argtypes[i] == ' ')
					r->argtypes[i] = b->argtypes[i];
				else if (b->argtypes[i] == ' ')
					r->argtypes[i] = a->argtypes[i];
				else if (a->argtypes[i] == b->argtypes[i])
					r->argtypes[i] = a->argtypes[i];
				else {
					PyErr_SetString(ErrorObject,
						"inconsistent argument lists");
					error = 1;
				}
			}
			for (; i < r->nargs; i++)
				if (i < a->nargs)
					r->argtypes[i] = a->argtypes[i];
				else
					r->argtypes[i] = b->argtypes[i];
			r->argtypes[i] = '\0';
		}
		else
			error = 1;
		if ((r->args = PyMem_NEW(stack_entry, r->nargs)) == NULL)
			error = 1;
		if ((r->constants = PyMem_NEW(stack_entry, r->nconstants))
				!= NULL) {
			for (i = 0; i < a->nconstants; i++)
				r->constants[i] = a->constants[i];
			for (j = 0; j < b->nconstants; i++, j++)
				r->constants[i] = b->constants[j];
		}
		else
			error = 1;
		if ((r->code = PyMem_NEW(int, r->codesize)) != NULL) {
			cp = r->code;
			for (i = 0; i < a->codesize; i++)
				*cp++ = a->code[i];
			if (a->type != type)
				*cp++ = TYPE_CAST + a->type + NTYPES*type;
			for (i = 0; i < b->codesize;) {
				*cp++ = j = b->code[i++];
				if (j == GET_CONSTANT) {
					*cp++ = b->code[i++] + a->nconstants;
				}
				else if (j == GET_ARG) {
					*cp++ = b->code[i++];
				}
			}
			if (b->type != type)
				*cp++ = TYPE_CAST + b->type + NTYPES*type;
			*cp = op;
		}
		else
			error = 1;
		if ((r->stack = PyMem_NEW(stack_entry, r->stacksize)) == NULL)
			error = 1;
		if (error) {
			numexpr_dealloc(r);
			r = NULL;
		}
	}
	return r;
}

static PyNumExprObject *
numexpr_unary_op(a, op)
	PyNumExprObject *a;
	int op;
{
	PyNumExprObject *r = new_numexprobject();
	int error = 0;
	int i;
	int *cp;

	if (r != NULL) {
		r->type = a->type;
		r->nargs = a->nargs;
		r->nconstants = a->nconstants;
		r->stacksize = a->stacksize;
		r->codesize = a->codesize + 1;
		if ((r->argtypes = PyMem_NEW(char, r->nargs+1)) != NULL)
			strcpy(r->argtypes, a->argtypes);
		else
			error = 1;
		if ((r->args = PyMem_NEW(stack_entry, r->nargs)) == NULL)
			error = 1;
		if ((r->constants = PyMem_NEW(stack_entry, r->nconstants))
				!= NULL)
			for (i = 0; i < a->nconstants; i++)
				r->constants[i] = a->constants[i];
		else
			error = 1;
		if ((r->code = PyMem_NEW(int, r->codesize)) != NULL) {
			cp = r->code;
			for (i = 0; i < a->codesize; i++)
				*cp++ = a->code[i];
			*cp = op;
		}
		if ((r->stack = PyMem_NEW(stack_entry, r->stacksize)) == NULL)
			error = 1;
		if (error) {
			numexpr_dealloc(r);
			r = NULL;
		}
	}
	return r;
}

static PyObject *
numexpr_add(a, b)
	PyNumExprObject *a;
	PyNumExprObject *b;
{
	char type = numexpr_commontype(a->type, b->type);
	return (PyObject *)numexpr_binary_op(a, b, type, ADD+type);
}

static PyObject *
numexpr_sub(a, b)
	PyNumExprObject *a;
	PyNumExprObject *b;
{
	char type = numexpr_commontype(a->type, b->type);
	return (PyObject *)numexpr_binary_op(a, b, type, SUB+type);
}

static PyObject *
numexpr_mul(a, b)
	PyNumExprObject *a;
	PyNumExprObject *b;
{
	char type = numexpr_commontype(a->type, b->type);
	return (PyObject *)numexpr_binary_op(a, b, type, MUL+type);
}

static PyObject *
numexpr_div(a, b)
	PyNumExprObject *a;
	PyNumExprObject *b;
{
	char type = numexpr_commontype(a->type, b->type);
	return (PyObject *)numexpr_binary_op(a, b, type, DIV+type);
}

static PyObject *
numexpr_pow(a, b)
	PyNumExprObject *a;
	PyNumExprObject *b;
{
	char type = numexpr_commontype(a->type, b->type);
	return (PyObject *)numexpr_binary_op(a, b, type, POW+type);
}

static PyObject *
numexpr_neg(a)
	PyNumExprObject *a;
{
	return (PyObject *)numexpr_unary_op(a, NEG+a->type);
}

static PyObject *
numexpr_pos(a)
	PyNumExprObject *a;
{
	Py_INCREF(a);
	return (PyObject *)a;
}

static PyObject *
numexpr_abs(a)
	PyNumExprObject *a;
{
	PyNumExprObject *r = numexpr_unary_op(a, ABS+a->type);
	if (r != NULL && r->type == COMPLEX)
		r->type = FLOAT;
	return (PyObject *)r;
}

static int
numexpr_nonzero(a)
	PyNumExprObject *a;
{
	return 1;
}


/* Type conversion */

static void
numexpr_convert(value, type_in, type_out)
	stack_entry *value;
	int type_in;
	int type_out;
{
	switch (type_in) {
		case INTEGER:
			if (type_out == FLOAT)
				value->f = value->i;
			else
				value->c.real = value->i;
				value->c.imag = 0;
			break;
		case FLOAT:
			if (type_out == INTEGER)
				value->i = value->f;
			else
				value->c.real = value->f;
				value->c.imag = 0;
			break;
		case COMPLEX:
			if (type_out == INTEGER)
				value->i = value->c.real;
			else
				value->f = value->c.real;
			break;
	}
}

/* Type coercion: Numbers are converted to constant expressions */

static PyObject *
numexpr_const(value, type_in, type_out)
	stack_entry value;
	int type_in, type_out;
{
	PyNumExprObject *rv;
	int error = 0;

	rv = new_numexprobject();
	if (rv != NULL) {
		rv->nargs = 0;
		rv->nconstants = 1;
		rv->codesize = 2;
		rv->stacksize = 1;
		rv->type = type_out;
		if (type_in != type_out)
			numexpr_convert(&value, type_in, type_out);
		if ((rv->constants = PyMem_NEW(stack_entry, rv->nconstants))
				!= NULL)
			*rv->constants = value;
		else
			error = 1;
		if ((rv->code = PyMem_NEW(int, rv->codesize)) != NULL) {
			rv->code[0] = GET_CONSTANT;
			rv->code[1] = 0;
		}
		else
			error = 1;
		if ((rv->stack = PyMem_NEW(stack_entry, rv->stacksize)) == NULL)
			error = 1;
		if (error) {
			numexpr_dealloc(rv);
			rv = NULL;
		}
	}
	return (PyObject *)rv;
}

static void
numexpr_coercenumber(value, type, number, numexpr)
	stack_entry value;
	int type;
	PyNumExprObject **number;
	PyNumExprObject **numexpr;
{
	int result_type = numexpr_commontype(type, (*numexpr)->type);

	*number = (PyNumExprObject *)numexpr_const(value, type, result_type);
	if ((*numexpr)->type != result_type)
		*numexpr = numexpr_unary_op(*numexpr,
			TYPE_CAST + (*numexpr)->type + NTYPES*result_type);
	else
		Py_INCREF(*numexpr);
}

static int
numexpr_coerce(a, b)
	PyObject **a;
	PyObject **b;
{
	PyObject **number;
	PyNumExprObject **numexpr;
	stack_entry value;

	if (PyNumExpr_Check(*a)) {
		numexpr = (PyNumExprObject **)a;
		number = b;
	}
	else {
		numexpr = (PyNumExprObject **)b;
		number = a;
	}
	if (PyInt_Check(*number)) {
		value.i = PyInt_AsLong(*number);
		numexpr_coercenumber(value, INTEGER, number, numexpr);
		return 0;
	}
	else if (PyFloat_Check(*number)) {
		value.f = PyFloat_AsDouble(*number);
		numexpr_coercenumber(value, FLOAT, number, numexpr);
		return 0;
	}
	else if (PyComplex_Check(*number)) {
		value.c = PyComplex_AsCComplex(*number);
		numexpr_coercenumber(value, COMPLEX, number, numexpr);
		return 0;
	}
	else
		return 1;
}

/* Module documentation string */

static char PyNumExpr_Type__doc__[] = 
"Compiled numerical expressions";

/* Type definition */

static PyNumberMethods numexpr_as_number = {
	(binaryfunc)numexpr_add, /*nb_add*/
	(binaryfunc)numexpr_sub, /*nb_subtract*/
	(binaryfunc)numexpr_mul, /*nb_multiply*/
	(binaryfunc)numexpr_div, /*nb_divide*/
	0,		/*nb_remainder*/
	0,		/*nb_divmod*/
	(ternaryfunc)numexpr_pow, /*nb_power*/
	(unaryfunc)numexpr_neg, /*nb_negative*/
	(unaryfunc)numexpr_pos, /*nb_positive*/
	(unaryfunc)numexpr_abs, /*nb_absolute*/
	(inquiry)numexpr_nonzero, /*nb_nonzero*/
	0,		/*nb_invert*/
	0,		/*nb_lshift*/
	0,		/*nb_rshift*/
	0,		/*nb_and*/
	0,		/*nb_xor*/
	0,		/*nb_or*/
	(coercion)numexpr_coerce, /*nb_coerce*/
	0,		/*nb_int*/
	0,		/*nb_long*/
	0,		/*nb_float*/
	0,		/*nb_oct*/
	0,		/*nb_hex*/
};

static PyTypeObject PyNumExpr_Type = {
	PyObject_HEAD_INIT(&PyType_Type)
	0,			/*ob_size*/
	"numexpr",			/*tp_name*/
	sizeof(PyNumExprObject),	/*tp_basicsize*/
	0,			/*tp_itemsize*/
	/* methods */
	(destructor)numexpr_dealloc, /*tp_dealloc*/
	0,			/*tp_print*/
	(getattrfunc)numexpr_getattr, /*tp_getattr*/
	0, 			/*tp_setattr*/
	0,			/*tp_compare*/
	(reprfunc)numexpr_repr, /*tp_repr*/
	&numexpr_as_number,	/*tp_as_number*/
	0,			/*tp_as_sequence*/
	0,			/*tp_as_mapping*/
	0,			/*tp_hash*/
	(ternaryfunc)numexpr_call,	/*tp_call*/
	(reprfunc)numexpr_repr,		/*tp_str*/

	/* Space for future expansion */
	0L,0L,0L,0L,
	PyNumExpr_Type__doc__	/* Documentation string */
};
/* --------------------------------------------------------------------- */

/* Return a new variable object. The (integer) argument
   indicates the position of the variable in the argument
   list of the compiled expression; it defaults to zero. */

static PyNumExprObject *
new_variable(args)
	PyObject *args;
{
	PyNumExprObject *rv;
	int argnum = 0;
	int error = 0;
	int i;

	PyArg_ParseTuple(args, "|i:numexpr variable", &argnum);
	rv = new_numexprobject();
	if (rv != NULL) {
		rv->nargs = argnum + 1;
		rv->stacksize = 1;
		rv->codesize = 2;
		if ((rv->argtypes = PyMem_NEW(char, rv->nargs+1)) != NULL) {
			for (i = 0; i < rv->nargs; i++)
				rv->argtypes[i] = ' ';
			rv->argtypes[i] = '\0';
		}
		else
			error = 1;
		if ((rv->args = PyMem_NEW(stack_entry, rv->nargs)) == NULL)
			error = 1;
		if ((rv->code = PyMem_NEW(int, rv->codesize)) != NULL) {
			rv->code[0] = GET_ARG;
			rv->code[1] = argnum;
		}
		else
			error = 1;
		if ((rv->stack = PyMem_NEW(stack_entry, rv->stacksize)) == NULL)
			error = 1;
		if (error) {
			numexpr_dealloc(rv);
			rv = NULL;
		}
	}
	return rv;
}

static PyObject *
numexpr_intvar(self, args)
	PyObject *self; /* Not used */
	PyObject *args;
{
	PyNumExprObject *rv = new_variable(args);
	if (rv != NULL) {
		rv->argtypes[rv->nargs-1] = 'l';
		rv->type = INTEGER;
	}
	return (PyObject *)rv;
}

static PyObject *
numexpr_floatvar(self, args)
	PyObject *self; /* Not used */
	PyObject *args;
{
	PyNumExprObject *rv = new_variable(args);
	if (rv != NULL) {
		rv->argtypes[rv->nargs-1] = 'd';
		rv->type = FLOAT;
	}
	return (PyObject *)rv;
}

static PyObject *
numexpr_complexvar(self, args)
	PyObject *self; /* Not used */
	PyObject *args;
{
	PyNumExprObject *rv = new_variable(args);
	if (rv != NULL) {
		rv->argtypes[rv->nargs-1] = 'D';
		rv->type = COMPLEX;
	}
	return (PyObject *)rv;
}

/* List of functions defined in the module */

static struct PyMethodDef numexpr_functions[] = {
	{"FloatVariable",	numexpr_floatvar, 1},
	{"IntegerVariable",	numexpr_intvar,	1},
	{"ComplexVariable",	numexpr_complexvar, 1},
	{NULL,	NULL}		/* sentinel */
};

/* Various math functions */

#if 0
static PyObject *
numexpr_sqrt(self)
	PyNumExprObject *self;
{
	PyNumExprObject *temp = self;
	PyNumExprObject *r;

	if (self->type == INTEGER) {
		temp = numexpr_unary_op(self, INTEGER_FLOAT);
		temp->type = FLOAT;
	}
	if (temp->type == COMPLEX)
		r = numexpr_unary_op(temp, SQRT_COMPLEX);
	else
		r = numexpr_unary_op(temp, SQRT_FLOAT);
	if (temp != self)
		numexpr_dealloc(temp);
	return (PyObject *)r;
}
#endif

static PyObject *
numexpr_mathfunc(arg, float_op, complex_op)
	PyNumExprObject *arg;
	int float_op;
	int complex_op;
{
	PyNumExprObject *temp = arg;
	PyNumExprObject *r;

	if (arg->type == INTEGER) {
		temp = numexpr_unary_op(arg, INTEGER_FLOAT);
		temp->type = FLOAT;
	}
	if (temp->type == COMPLEX)
		r = numexpr_unary_op(temp, complex_op);
	else
		r = numexpr_unary_op(temp, float_op);
	if (temp != arg)
		numexpr_dealloc(temp);
	return (PyObject *)r;
}

#define add_mathfunc(func_name, float_op, complex_op) \
	static PyObject * \
	func_name(self) \
		PyNumExprObject *self; \
	{ \
		return numexpr_mathfunc(self, float_op, complex_op); \
	}

add_mathfunc(numexpr_acos, ACOS_FLOAT, ACOS_COMPLEX)
add_mathfunc(numexpr_acosh, ACOSH_FLOAT, ACOSH_COMPLEX)
add_mathfunc(numexpr_asin, ASIN_FLOAT, ASIN_COMPLEX)
add_mathfunc(numexpr_asinh, ASINH_FLOAT, ASINH_COMPLEX)
add_mathfunc(numexpr_atan, ATAN_FLOAT, ATAN_COMPLEX)
add_mathfunc(numexpr_atanh, ATANH_FLOAT, ATANH_COMPLEX)
add_mathfunc(numexpr_cos, COS_FLOAT, COS_COMPLEX)
add_mathfunc(numexpr_cosh, COSH_FLOAT, COSH_COMPLEX)
add_mathfunc(numexpr_exp, EXP_FLOAT, EXP_COMPLEX)
add_mathfunc(numexpr_log, LOG_FLOAT, LOG_COMPLEX)
add_mathfunc(numexpr_log10, LOG10_FLOAT, LOG10_COMPLEX)
add_mathfunc(numexpr_sin, SIN_FLOAT, SIN_COMPLEX)
add_mathfunc(numexpr_sinh, SINH_FLOAT, SINH_COMPLEX)
add_mathfunc(numexpr_sqrt, SQRT_FLOAT, SQRT_COMPLEX)
add_mathfunc(numexpr_tan, TAN_FLOAT, TAN_COMPLEX)
add_mathfunc(numexpr_tanh, TANH_FLOAT, TANH_COMPLEX)


/* List of methods of type "numexpr" */

static struct PyMethodDef numexpr_methods[] = {
	{"acos",	(PyCFunction)numexpr_acos, 1},
	{"acosh",	(PyCFunction)numexpr_acosh, 1},
	{"asin",	(PyCFunction)numexpr_asin, 1},
	{"asinh",	(PyCFunction)numexpr_asinh, 1},
	{"atan",	(PyCFunction)numexpr_atan, 1},
	{"atanh",	(PyCFunction)numexpr_atanh, 1},
	{"cos",		(PyCFunction)numexpr_cos, 1},
	{"cosh",	(PyCFunction)numexpr_cosh, 1},
	{"exp",		(PyCFunction)numexpr_exp, 1},
	{"log",		(PyCFunction)numexpr_log, 1},
	{"log10",	(PyCFunction)numexpr_log10, 1},
	{"sin",		(PyCFunction)numexpr_sin, 1},
	{"sinh",	(PyCFunction)numexpr_sinh, 1},
	{"sqrt",	(PyCFunction)numexpr_sqrt, 1},
	{"tan",		(PyCFunction)numexpr_tan, 1},
	{"tanh",	(PyCFunction)numexpr_tanh, 1},
	{NULL,		NULL}		/* sentinel */
};


/* Initialization function for the module (*must* be called initnumexpr) */

void
initnumexpr()
{
	PyObject *m, *d;

	/* Create the module and add the functions */
	m = Py_InitModule("numexpr", numexpr_functions);

	/* Add some symbolic constants to the module */
	d = PyModule_GetDict(m);
	ErrorObject = PyString_FromString("NumExprError");
	PyDict_SetItemString(d, "error", ErrorObject);

	/* Check for errors */
	if (PyErr_Occurred())
		Py_FatalError("can't initialize module numexpr");
}

=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================