Daniel@0: /*
Daniel@0: 
Daniel@0: 	GENOPS.C
Daniel@0: 		Generalized arithmetic operators overloading built-in functions.
Daniel@0: 
Daniel@0: 	written by Douglas M. Schwarz
Daniel@0: 	schwarz@servtech.com
Daniel@0: 	26 December 1998
Daniel@0: 	Last modified: 2 April 1999
Daniel@0: 
Daniel@0: 	Copyright 1998-1999 by Douglas M. Schwarz.  All rights reserved.
Daniel@0: 
Daniel@0: */
Daniel@0: 
Daniel@0: 
Daniel@0: /*
Daniel@0: 
Daniel@0: Build MEX file by entering the appropriate command at the MATLAB prompt
Daniel@0: (-D<name> option is equivalent to #define <name> in source file):
Daniel@0: 
Daniel@0: mex genops.c -DPLUS_MEX -output plus
Daniel@0: mex genops.c -DMINUS_MEX -output minus
Daniel@0: mex genops.c -DTIMES_MEX -output times
Daniel@0: mex genops.c -DRDIVIDE_MEX -output rdivide
Daniel@0: mex genops.c -DLDIVIDE_MEX -output ldivide
Daniel@0: mex genops.c -DPOWER_MEX -output power
Daniel@0: mex genops.c -DEQ_MEX -output eq
Daniel@0: mex genops.c -DNE_MEX -output ne
Daniel@0: mex genops.c -DLT_MEX -output lt
Daniel@0: mex genops.c -DGT_MEX -output gt
Daniel@0: mex genops.c -DLE_MEX -output le
Daniel@0: mex genops.c -DGE_MEX -output ge
Daniel@0: 
Daniel@0: */
Daniel@0: 
Daniel@0: /*	This file has been formatted for a tab equal to 4 spaces. */
Daniel@0: 
Daniel@0: #if defined(EQ_MEX) || defined(NE_MEX) || defined(LT_MEX) || defined(GT_MEX) \
Daniel@0: 		|| defined(LE_MEX) || defined(GE_MEX)
Daniel@0: #define	RELOP_MEX
Daniel@0: #endif
Daniel@0: 
Daniel@0: #include "mex.h"
Daniel@0: #include "matrix.h"
Daniel@0: #ifdef POWER_MEX
Daniel@0: #include <math.h>
Daniel@0: #define PI 3.141592653589793
Daniel@0: #endif
Daniel@0: 
Daniel@0: bool allequal(int, const int *, const int *);
Daniel@0: void removeZeroImag(double *, double *, int, const int *, int, mxArray **);
Daniel@0: 
Daniel@0: #define	xMat  prhs[0]
Daniel@0: #define	yMat  prhs[1]
Daniel@0: #define	zMat  plhs[0]
Daniel@0: 
Daniel@0: #define	min(A,B)  ((A) < (B) ? (A) : (B))
Daniel@0: #define	max(A,B)  ((A) > (B) ? (A) : (B))
Daniel@0: 
Daniel@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
Daniel@0: {
Daniel@0: 	double		*xrp, *xip, *yrp, *yip;
Daniel@0: #ifndef RELOP_MEX
Daniel@0: 	double		*zr, *zi, *zip;
Daniel@0: #endif
Daniel@0: 	double		*zrp, *zrend;
Daniel@0: 	int			xnd, ynd, numElements = 1;
Daniel@0: 	const int	*xdim, *ydim;
Daniel@0: 	bool		xcmplx, ycmplx;
Daniel@0: 	mxClassID	yclass;
Daniel@0: 	int			*s, ndim, *sx, *sy, i, *cpsx, *cpsy;
Daniel@0: 	int			*subs, *s1, *cpsx2, *cpsy2;
Daniel@0: 	int			ix = 0, iy = 0;
Daniel@0: 	mxArray		*args[3], *result[1];
Daniel@0: #if defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
Daniel@0: 	double		denom;
Daniel@0: #endif
Daniel@0: #ifdef POWER_MEX
Daniel@0: 	double		mag, theta, phi, magx;
Daniel@0: 	int			flops = 0;
Daniel@0: #endif
Daniel@0: 	
Daniel@0: 	
Daniel@0: 	if (nrhs != 2)
Daniel@0: 		mexErrMsgTxt("Incorrect number of inputs.");
Daniel@0: 	
Daniel@0: 	if (nlhs > 1)
Daniel@0: 		mexErrMsgTxt("Too many output arguments.");
Daniel@0: 	
Daniel@0: 	xnd = mxGetNumberOfDimensions(xMat);
Daniel@0: 	ynd = mxGetNumberOfDimensions(yMat);
Daniel@0: 	xdim = mxGetDimensions(xMat);
Daniel@0: 	ydim = mxGetDimensions(yMat);
Daniel@0: 	
Daniel@0: 	yclass = mxGetClassID(yMat);
Daniel@0: 	
Daniel@0: /*	If the built-in function in MATLAB can handle the arguments
Daniel@0: 	then use that. */
Daniel@0: 	if (yclass != mxDOUBLE_CLASS || 
Daniel@0: 		(xnd == 2  &&  xdim[0] == 1  &&  xdim[1] == 1) || 
Daniel@0: 		(ynd == 2  &&  ydim[0] == 1  &&  ydim[1] == 1) || 
Daniel@0: 		(xnd == ynd  &&  allequal(xnd,xdim,ydim)))
Daniel@0: 	{
Daniel@0: #ifdef PLUS_MEX
Daniel@0: 		args[0] = mxCreateString("plus");
Daniel@0: #elif defined(MINUS_MEX)
Daniel@0: 		args[0] = mxCreateString("minus");
Daniel@0: #elif defined(TIMES_MEX)
Daniel@0: 		args[0] = mxCreateString("times");
Daniel@0: #elif defined(RDIVIDE_MEX)
Daniel@0: 		args[0] = mxCreateString("rdivide");
Daniel@0: #elif defined(LDIVIDE_MEX)
Daniel@0: 		args[0] = mxCreateString("ldivide");
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 		args[0] = mxCreateString("power");
Daniel@0: #elif defined(EQ_MEX)
Daniel@0: 		args[0] = mxCreateString("eq");
Daniel@0: #elif defined(NE_MEX)
Daniel@0: 		args[0] = mxCreateString("ne");
Daniel@0: #elif defined(LT_MEX)
Daniel@0: 		args[0] = mxCreateString("lt");
Daniel@0: #elif defined(GT_MEX)
Daniel@0: 		args[0] = mxCreateString("gt");
Daniel@0: #elif defined(LE_MEX)
Daniel@0: 		args[0] = mxCreateString("le");
Daniel@0: #elif defined(GE_MEX)
Daniel@0: 		args[0] = mxCreateString("ge");
Daniel@0: #endif
Daniel@0: 		args[1] = (mxArray *)xMat;
Daniel@0: 		args[2] = (mxArray *)yMat;
Daniel@0: 		mexCallMATLAB(1, result, 3, args, "builtin");
Daniel@0: 		mxDestroyArray(args[0]);
Daniel@0: 		zMat = result[0];
Daniel@0: 	}
Daniel@0: 	else  /* X and Y are both N-D and different dimensionality. */
Daniel@0: 	{
Daniel@0: 		ndim = max(xnd,ynd);
Daniel@0: 		sx = (int *)mxMalloc(sizeof(int)*ndim);
Daniel@0: 		sy = (int *)mxMalloc(sizeof(int)*ndim);
Daniel@0: 		s =  (int *)mxMalloc(sizeof(int)*ndim);
Daniel@0: 		s1 = (int *)mxMalloc(sizeof(int)*ndim);
Daniel@0: 		*(cpsx = (int *)mxMalloc(sizeof(int)*ndim)) = 1;
Daniel@0: 		*(cpsy = (int *)mxMalloc(sizeof(int)*ndim)) = 1;
Daniel@0: 		subs =   (int *)mxMalloc(sizeof(int)*ndim);
Daniel@0: 		cpsx2 =  (int *)mxMalloc(sizeof(int)*ndim);
Daniel@0: 		cpsy2 =  (int *)mxMalloc(sizeof(int)*ndim);
Daniel@0: 		for (i = 0; i < ndim; i++)
Daniel@0: 		{
Daniel@0: 			subs[i] = 0;
Daniel@0: 			sx[i] = (i < xnd) ? xdim[i] : 1;
Daniel@0: 			sy[i] = (i < ynd) ? ydim[i] : 1;
Daniel@0: 			if (sx[i] == sy[i])
Daniel@0: 				s[i] = sx[i];
Daniel@0: 			else if (sx[i] == 1)
Daniel@0: 				s[i] = sy[i];
Daniel@0: 			else if (sy[i] == 1)
Daniel@0: 				s[i] = sx[i];
Daniel@0: 			else
Daniel@0: 			{
Daniel@0: 				mxFree(sx);
Daniel@0: 				mxFree(sy);
Daniel@0: 				mxFree(s);
Daniel@0: 				mxFree(s1);
Daniel@0: 				mxFree(cpsx);
Daniel@0: 				mxFree(cpsy);
Daniel@0: 				mxFree(subs);
Daniel@0: 				mxFree(cpsx2);
Daniel@0: 				mxFree(cpsy2);
Daniel@0: 				mexErrMsgTxt("Array dimensions are not appropriate.");
Daniel@0: 			}
Daniel@0: 			s1[i] = s[i] - 1;
Daniel@0: 			numElements *= s[i];
Daniel@0: 		}
Daniel@0: 				
Daniel@0: 		for (i = 0; i < ndim-1; i++)
Daniel@0: 		{
Daniel@0: 			cpsx[i+1] = cpsx[i]*sx[i]--;
Daniel@0: 			cpsy[i+1] = cpsy[i]*sy[i]--;
Daniel@0: 			cpsx2[i] = cpsx[i]*sx[i];
Daniel@0: 			cpsy2[i] = cpsy[i]*sy[i];
Daniel@0: 		}
Daniel@0: 		cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
Daniel@0: 		cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
Daniel@0: 		
Daniel@0: 		xcmplx = mxIsComplex(xMat);
Daniel@0: 		ycmplx = mxIsComplex(yMat);
Daniel@0: 		
Daniel@0: 		if (!xcmplx && !ycmplx)  /* X and Y both N-D, both real. */
Daniel@0: 		{
Daniel@0: #ifdef POWER_MEX
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
Daniel@0: 			zrp = zr = mxGetPr(zMat);
Daniel@0: 			zip = zi = mxGetPi(zMat);
Daniel@0: #elif defined(RELOP_MEX)
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
Daniel@0: 			mxSetLogical(zMat);
Daniel@0: 			zrp = mxGetPr(zMat);
Daniel@0: #else
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
Daniel@0: 			zrp = mxGetPr(zMat);
Daniel@0: #endif
Daniel@0: 			xrp = mxGetPr(xMat);
Daniel@0: 			yrp = mxGetPr(yMat);
Daniel@0: 			zrend = zrp + numElements;
Daniel@0: 			while (zrp < zrend)
Daniel@0: 			{
Daniel@0: #ifdef PLUS_MEX
Daniel@0: 				*zrp++ = *xrp + *yrp;
Daniel@0: #elif defined(MINUS_MEX)
Daniel@0: 				*zrp++ = *xrp - *yrp;
Daniel@0: #elif defined(TIMES_MEX)
Daniel@0: 				*zrp++ = *xrp * *yrp;
Daniel@0: #elif defined(RDIVIDE_MEX)
Daniel@0: 				*zrp++ = *xrp / *yrp;
Daniel@0: #elif defined(LDIVIDE_MEX)
Daniel@0: 				*zrp++ = *yrp / *xrp;
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 				if (*xrp < 0.0 && *yrp != floor(*yrp))
Daniel@0: 				{
Daniel@0: 					mag = pow(-*xrp,*yrp);
Daniel@0: 					theta = PI * *yrp;
Daniel@0: 					*zrp++ = mag*cos(theta);
Daniel@0: 					*zip++ = mag*sin(theta);
Daniel@0: 					flops += 18;
Daniel@0: 				}
Daniel@0: 				else
Daniel@0: 				{
Daniel@0: 					*zrp++ = pow(*xrp,*yrp);
Daniel@0: 					*zip++ = 0.0;
Daniel@0: 					flops++;
Daniel@0: 				}
Daniel@0: #elif defined(EQ_MEX)
Daniel@0: 				*zrp++ = (*xrp == *yrp);
Daniel@0: #elif defined(NE_MEX)
Daniel@0: 				*zrp++ = (*xrp != *yrp);
Daniel@0: #elif defined(LT_MEX)
Daniel@0: 				*zrp++ = (*xrp < *yrp);
Daniel@0: #elif defined(GT_MEX)
Daniel@0: 				*zrp++ = (*xrp > *yrp);
Daniel@0: #elif defined(LE_MEX)
Daniel@0: 				*zrp++ = (*xrp <= *yrp);
Daniel@0: #elif defined(GE_MEX)
Daniel@0: 				*zrp++ = (*xrp >= *yrp);
Daniel@0: #endif
Daniel@0: 				for (i = 0; i < ndim; i++)
Daniel@0: 				{
Daniel@0: 					if (subs[i] == s1[i])
Daniel@0: 					{
Daniel@0: 						subs[i] = 0;
Daniel@0: 						if (sx[i])
Daniel@0: 							xrp -= cpsx2[i];
Daniel@0: 						if (sy[i])
Daniel@0: 							yrp -= cpsy2[i];
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						subs[i]++;
Daniel@0: 						if (sx[i])
Daniel@0: 							xrp += cpsx[i];
Daniel@0: 						if (sy[i])
Daniel@0: 							yrp += cpsy[i];
Daniel@0: 						break;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 			}
Daniel@0: #ifdef POWER_MEX
Daniel@0: 			mexAddFlops(flops);
Daniel@0: 			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
Daniel@0: #elif !defined(RELOP_MEX)
Daniel@0: 			mexAddFlops(numElements);
Daniel@0: #endif
Daniel@0: 		}
Daniel@0: 		else if (!xcmplx && ycmplx)  /* X and Y both N-D, X real, Y complex. */
Daniel@0: 		{
Daniel@0: #ifdef POWER_MEX
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
Daniel@0: 			zrp = zr = mxGetPr(zMat);
Daniel@0: 			zip = zi = mxGetPi(zMat);
Daniel@0: #elif defined(RELOP_MEX)
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
Daniel@0: 			mxSetLogical(zMat);
Daniel@0: 			zrp = mxGetPr(zMat);
Daniel@0: #else
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
Daniel@0: 			zrp = mxGetPr(zMat);
Daniel@0: 			zip = mxGetPi(zMat);
Daniel@0: #endif
Daniel@0: 			xrp = mxGetPr(xMat);
Daniel@0: 			yrp = mxGetPr(yMat);
Daniel@0: 			yip = mxGetPi(yMat);
Daniel@0: 			zrend = zrp + numElements;
Daniel@0: 			while (zrp < zrend)
Daniel@0: 			{
Daniel@0: #ifdef PLUS_MEX
Daniel@0: 				*zrp++ = *xrp + *yrp;
Daniel@0: 				*zip++ = *yip;
Daniel@0: #elif defined(MINUS_MEX)
Daniel@0: 				*zrp++ = *xrp - *yrp;
Daniel@0: 				*zip++ = -*yip;
Daniel@0: #elif defined(TIMES_MEX)
Daniel@0: 				*zrp++ = *xrp * *yrp;
Daniel@0: 				*zip++ = *xrp * *yip;
Daniel@0: #elif defined(RDIVIDE_MEX)
Daniel@0: 				denom = *yrp * *yrp + *yip * *yip;
Daniel@0: 				*zrp++ = (*xrp * *yrp)/denom;
Daniel@0: 				*zip++ = (-*xrp * *yip)/denom;
Daniel@0: #elif defined(LDIVIDE_MEX)
Daniel@0: 				*zrp++ = *yrp / *xrp;
Daniel@0: 				*zip++ = *yip / *xrp;
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 				if (*yip == 0.0)
Daniel@0: 				{
Daniel@0: 					if (*xrp < 0.0 && *yrp != floor(*yrp))
Daniel@0: 					{
Daniel@0: 						mag = pow(-*xrp,*yrp);
Daniel@0: 						theta = PI * *yrp;
Daniel@0: 						*zrp++ = mag*cos(theta);
Daniel@0: 						*zip++ = mag*sin(theta);
Daniel@0: 						flops += 18;
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						*zrp++ = pow(*xrp,*yrp);
Daniel@0: 						*zip++ = 0.0;
Daniel@0: 						flops++;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 				else
Daniel@0: 				{
Daniel@0: 					if (*xrp < 0.0)
Daniel@0: 					{
Daniel@0: 						mag = pow(-*xrp,*yrp)*exp(-PI * *yip);
Daniel@0: 						theta = *yip * log(-*xrp) + PI * *yrp;
Daniel@0: 						*zrp++ = mag*cos(theta);
Daniel@0: 						*zip++ = mag*sin(theta);
Daniel@0: 						flops += 18;
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						mag = pow(*xrp,*yrp);
Daniel@0: 						theta = *yip * log(*xrp);
Daniel@0: 						*zrp++ = mag*cos(theta);
Daniel@0: 						*zip++ = mag*sin(theta);
Daniel@0: 						flops += 13;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: #elif defined(EQ_MEX)
Daniel@0: 				*zrp++ = (*xrp == *yrp) && (*yip == 0.0);
Daniel@0: #elif defined(NE_MEX)
Daniel@0: 				*zrp++ = (*xrp != *yrp) || (*yip != 0.0);
Daniel@0: #elif defined(LT_MEX)
Daniel@0: 				*zrp++ = (*xrp < *yrp);
Daniel@0: #elif defined(GT_MEX)
Daniel@0: 				*zrp++ = (*xrp > *yrp);
Daniel@0: #elif defined(LE_MEX)
Daniel@0: 				*zrp++ = (*xrp <= *yrp);
Daniel@0: #elif defined(GE_MEX)
Daniel@0: 				*zrp++ = (*xrp >= *yrp);
Daniel@0: #endif
Daniel@0: 				for (i = 0; i < ndim; i++)
Daniel@0: 				{
Daniel@0: 					if (subs[i] == s1[i])
Daniel@0: 					{
Daniel@0: 						subs[i] = 0;
Daniel@0: 						if (sx[i])
Daniel@0: 							xrp -= cpsx2[i];
Daniel@0: 						if (sy[i])
Daniel@0: 						{
Daniel@0: 							yrp -= cpsy2[i];
Daniel@0: 							yip -= cpsy2[i];
Daniel@0: 						}
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						subs[i]++;
Daniel@0: 						if (sx[i])
Daniel@0: 							xrp += cpsx[i];
Daniel@0: 						if (sy[i])
Daniel@0: 						{
Daniel@0: 							yrp += cpsy[i];
Daniel@0: 							yip += cpsy[i];
Daniel@0: 						}
Daniel@0: 						break;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 			}
Daniel@0: #if defined(PLUS_MEX) || defined(MINUS_MEX)
Daniel@0: 			mexAddFlops(2*numElements);
Daniel@0: #elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
Daniel@0: 			mexAddFlops(6*numElements);
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 			mexAddFlops(flops);
Daniel@0: 			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
Daniel@0: #endif
Daniel@0: 		}
Daniel@0: 		else if (xcmplx && !ycmplx)  /* X and Y both N-D, X complex, Y real. */
Daniel@0: 		{
Daniel@0: #ifdef POWER_MEX
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
Daniel@0: 			zrp = zr = mxGetPr(zMat);
Daniel@0: 			zip = zi = mxGetPi(zMat);
Daniel@0: #elif defined(RELOP_MEX)
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
Daniel@0: 			mxSetLogical(zMat);
Daniel@0: 			zrp = mxGetPr(zMat);
Daniel@0: #else
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
Daniel@0: 			zrp = mxGetPr(zMat);
Daniel@0: 			zip = mxGetPi(zMat);
Daniel@0: #endif
Daniel@0: 			xrp = mxGetPr(xMat);
Daniel@0: 			xip = mxGetPi(xMat);
Daniel@0: 			yrp = mxGetPr(yMat);
Daniel@0: 			zrend = zrp + numElements;
Daniel@0: 			while (zrp < zrend)
Daniel@0: 			{
Daniel@0: #ifdef PLUS_MEX
Daniel@0: 				*zrp++ = *xrp + *yrp;
Daniel@0: 				*zip++ = *xip;
Daniel@0: #elif defined(MINUS_MEX)
Daniel@0: 				*zrp++ = *xrp - *yrp;
Daniel@0: 				*zip++ = *xip;
Daniel@0: #elif defined(TIMES_MEX)
Daniel@0: 				*zrp++ = *xrp * *yrp;
Daniel@0: 				*zip++ = *xip * *yrp;
Daniel@0: #elif defined(RDIVIDE_MEX)
Daniel@0: 				*zrp++ = *xrp / *yrp;
Daniel@0: 				*zip++ = *xip / *yrp;
Daniel@0: #elif defined(LDIVIDE_MEX)
Daniel@0: 				denom = *xrp * *xrp + *xip * *xip;
Daniel@0: 				*zrp++ = (*xrp * *yrp)/denom;
Daniel@0: 				*zip++ = (-*xip * *yrp)/denom;
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 				if (*xip == 0.0)
Daniel@0: 				{
Daniel@0: 					if (*xrp < 0.0 && *yrp != floor(*yrp))
Daniel@0: 					{
Daniel@0: 						mag = pow(-*xrp,*yrp);
Daniel@0: 						theta = PI * *yrp;
Daniel@0: 						*zrp++ = mag*cos(theta);
Daniel@0: 						*zip++ = mag*sin(theta);
Daniel@0: 						flops += 18;
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						*zrp++ = pow(*xrp,*yrp);
Daniel@0: 						*zip++ = 0.0;
Daniel@0: 						flops++;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 				else
Daniel@0: 				{
Daniel@0: 					mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp);
Daniel@0: 					theta = *yrp*atan2(*xip,*xrp);
Daniel@0: 					*zrp++ = mag*cos(theta);
Daniel@0: 					*zip++ = mag*sin(theta);
Daniel@0: 					flops += 18;
Daniel@0: 				}
Daniel@0: #elif defined(EQ_MEX)
Daniel@0: 				*zrp++ = (*xrp == *yrp) && (*xip == 0.0);
Daniel@0: #elif defined(NE_MEX)
Daniel@0: 				*zrp++ = (*xrp != *yrp) || (*xip != 0.0);
Daniel@0: #elif defined(LT_MEX)
Daniel@0: 				*zrp++ = (*xrp < *yrp);
Daniel@0: #elif defined(GT_MEX)
Daniel@0: 				*zrp++ = (*xrp > *yrp);
Daniel@0: #elif defined(LE_MEX)
Daniel@0: 				*zrp++ = (*xrp <= *yrp);
Daniel@0: #elif defined(GE_MEX)
Daniel@0: 				*zrp++ = (*xrp >= *yrp);
Daniel@0: #endif
Daniel@0: 				for (i = 0; i < ndim; i++)
Daniel@0: 				{
Daniel@0: 					if (subs[i] == s1[i])
Daniel@0: 					{
Daniel@0: 						subs[i] = 0;
Daniel@0: 						if (sx[i])
Daniel@0: 						{
Daniel@0: 							xrp -= cpsx2[i];
Daniel@0: 							xip -= cpsx2[i];
Daniel@0: 						}
Daniel@0: 						if (sy[i])
Daniel@0: 							yrp -= cpsy2[i];
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						subs[i]++;
Daniel@0: 						if (sx[i])
Daniel@0: 						{
Daniel@0: 							xrp += cpsx[i];
Daniel@0: 							xip += cpsx[i];
Daniel@0: 						}
Daniel@0: 						if (sy[i])
Daniel@0: 							yrp += cpsy[i];
Daniel@0: 						break;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 			}
Daniel@0: #if defined(PLUS_MEX) || defined(MINUS_MEX)
Daniel@0: 			mexAddFlops(2*numElements);
Daniel@0: #elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
Daniel@0: 			mexAddFlops(6*numElements);
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 			mexAddFlops(flops);
Daniel@0: 			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
Daniel@0: #endif
Daniel@0: 		}
Daniel@0: 		else if (xcmplx && ycmplx)  /* X and Y both N-D, both complex. */
Daniel@0: 		{
Daniel@0: #if defined(RELOP_MEX)
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
Daniel@0: 			mxSetLogical(zMat);
Daniel@0: 			zrp = mxGetPr(zMat);
Daniel@0: #else
Daniel@0: 			zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
Daniel@0: 			zrp = zr = mxGetPr(zMat);
Daniel@0: 			zip = zi = mxGetPi(zMat);
Daniel@0: #endif
Daniel@0: 			xrp = mxGetPr(xMat);
Daniel@0: 			xip = mxGetPi(xMat);
Daniel@0: 			yrp = mxGetPr(yMat);
Daniel@0: 			yip = mxGetPi(yMat);
Daniel@0: 			zrend = zrp + numElements;
Daniel@0: 			while (zrp < zrend)
Daniel@0: 			{
Daniel@0: #ifdef PLUS_MEX
Daniel@0: 				*zrp++ = *xrp + *yrp;
Daniel@0: 				*zip++ = *xip + *yip;
Daniel@0: #elif defined(MINUS_MEX)
Daniel@0: 				*zrp++ = *xrp - *yrp;
Daniel@0: 				*zip++ = *xip - *yip;
Daniel@0: #elif defined(TIMES_MEX)
Daniel@0: 				*zrp++ = *xrp * *yrp - *xip * *yip;
Daniel@0: 				*zip++ = *xip * *yrp + *xrp * *yip;
Daniel@0: #elif defined(RDIVIDE_MEX)
Daniel@0: 				denom = *yrp * *yrp + *yip * *yip;
Daniel@0: 				*zrp++ = (*xrp * *yrp + *xip * *yip)/denom;
Daniel@0: 				*zip++ = (*xip * *yrp - *xrp * *yip)/denom;
Daniel@0: #elif defined(LDIVIDE_MEX)
Daniel@0: 				denom = *xrp * *xrp + *xip * *xip;
Daniel@0: 				*zrp++ = (*xrp * *yrp + *xip * *yip)/denom;
Daniel@0: 				*zip++ = (*xrp * *yip - *xip * *yrp)/denom;
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 				if (*xip == 0.0 && *yip == 0.0)
Daniel@0: 				{
Daniel@0: 					if (*xrp < 0.0 && *yrp != floor(*yrp))
Daniel@0: 					{
Daniel@0: 						mag = pow(-*xrp,*yrp);
Daniel@0: 						theta = PI * *yrp;
Daniel@0: 						*zrp++ = mag*cos(theta);
Daniel@0: 						*zip++ = mag*sin(theta);
Daniel@0: 						flops += 18;
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						*zrp++ = pow(*xrp,*yrp);
Daniel@0: 						*zip++ = 0.0;
Daniel@0: 						flops++;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 				else if (*xip == 0.0)
Daniel@0: 				{
Daniel@0: 					if (*xrp < 0.0)
Daniel@0: 					{
Daniel@0: 						mag = pow(-*xrp,*yrp)*exp(-PI * *yip);
Daniel@0: 						theta = *yip * log(-*xrp) + PI * *yrp;
Daniel@0: 						*zrp++ = mag*cos(theta);
Daniel@0: 						*zip++ = mag*sin(theta);
Daniel@0: 						flops += 18;
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						mag = pow(*xrp,*yrp);
Daniel@0: 						theta = *yip * log(*xrp);
Daniel@0: 						*zrp++ = mag*cos(theta);
Daniel@0: 						*zip++ = mag*sin(theta);
Daniel@0: 						flops += 13;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 				else if (*yip == 0.0)
Daniel@0: 				{
Daniel@0: 					mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp);
Daniel@0: 					theta = *yrp * atan2(*xip,*xrp);
Daniel@0: 					*zrp++ = mag*cos(theta);
Daniel@0: 					*zip++ = mag*sin(theta);
Daniel@0: 					flops += 18;
Daniel@0: 				}
Daniel@0: 				else
Daniel@0: 				{
Daniel@0: 					magx = sqrt(*xrp * *xrp + *xip * *xip);
Daniel@0: 					phi = atan2(*xip,*xrp);
Daniel@0: 					mag = pow(magx,*yrp)*exp(-*yip * phi);
Daniel@0: 					theta = *yip * log(magx) + *yrp * phi;
Daniel@0: 					*zrp++ = mag*cos(theta);
Daniel@0: 					*zip++ = mag*sin(theta);
Daniel@0: 					flops += 18;
Daniel@0: 				}
Daniel@0: #elif defined(EQ_MEX)
Daniel@0: 				*zrp++ = (*xrp == *yrp) && (*xip == *yip);
Daniel@0: #elif defined(NE_MEX)
Daniel@0: 				*zrp++ = (*xrp != *yrp) || (*xip != *yip);
Daniel@0: #elif defined(LT_MEX)
Daniel@0: 				*zrp++ = (*xrp < *yrp);
Daniel@0: #elif defined(GT_MEX)
Daniel@0: 				*zrp++ = (*xrp > *yrp);
Daniel@0: #elif defined(LE_MEX)
Daniel@0: 				*zrp++ = (*xrp <= *yrp);
Daniel@0: #elif defined(GE_MEX)
Daniel@0: 				*zrp++ = (*xrp >= *yrp);
Daniel@0: #endif
Daniel@0: 				for (i = 0; i < ndim; i++)
Daniel@0: 				{
Daniel@0: 					if (subs[i] == s1[i])
Daniel@0: 					{
Daniel@0: 						subs[i] = 0;
Daniel@0: 						if (sx[i])
Daniel@0: 						{
Daniel@0: 							xrp -= cpsx2[i];
Daniel@0: 							xip -= cpsx2[i];
Daniel@0: 						}
Daniel@0: 						if (sy[i])
Daniel@0: 						{
Daniel@0: 							yrp -= cpsy2[i];
Daniel@0: 							yip -= cpsy2[i];
Daniel@0: 						}
Daniel@0: 					}
Daniel@0: 					else
Daniel@0: 					{
Daniel@0: 						subs[i]++;
Daniel@0: 						if (sx[i])
Daniel@0: 						{
Daniel@0: 							xrp += cpsx[i];
Daniel@0: 							xip += cpsx[i];
Daniel@0: 						}
Daniel@0: 						if (sy[i])
Daniel@0: 						{
Daniel@0: 							yrp += cpsy[i];
Daniel@0: 							yip += cpsy[i];
Daniel@0: 						}
Daniel@0: 						break;
Daniel@0: 					}
Daniel@0: 				}
Daniel@0: 			}
Daniel@0: #if defined(PLUS_MEX) || defined(MINUS_MEX)
Daniel@0: 			mexAddFlops(2*numElements);
Daniel@0: #elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
Daniel@0: 			mexAddFlops(6*numElements);
Daniel@0: #elif defined(POWER_MEX)
Daniel@0: 			mexAddFlops(flops);
Daniel@0: #endif
Daniel@0: #ifndef RELOP_MEX
Daniel@0: 			removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
Daniel@0: #endif
Daniel@0: 		}
Daniel@0: 	}
Daniel@0: }
Daniel@0: 
Daniel@0: 
Daniel@0: /***********************************************************
Daniel@0: *                                                          *
Daniel@0: *   Tests to see if the vectors xdim and ydim are equal.   *
Daniel@0: *                                                          *
Daniel@0: ***********************************************************/
Daniel@0: bool allequal(int ndim, const int *xdim, const int *ydim)
Daniel@0: {
Daniel@0: 	int		i;
Daniel@0: 	bool	result = true;
Daniel@0: 	
Daniel@0: 	for (i = 0; i < ndim; i++)
Daniel@0: 		result = result && (xdim[i] == ydim[i]);
Daniel@0: 	
Daniel@0: 	return(result);
Daniel@0: }
Daniel@0: 
Daniel@0: 
Daniel@0: /******************************************************************************
Daniel@0: *                                                                             *
Daniel@0: *   Tests to see if every imaginary element is identically zero and, if so,   *
Daniel@0: *   creates a new array which is real and copies the real elements to it.     *
Daniel@0: *                                                                             *
Daniel@0: ******************************************************************************/
Daniel@0: void removeZeroImag(double *zr, double *zi, int ndim, const int *s,
Daniel@0: 					int numElements, mxArray *plhs[])
Daniel@0: {
Daniel@0: 	double			*zrend, *ziend, *zip, *z1p, *z2p;
Daniel@0: 	bool			allImZero = true;
Daniel@0: 	mxArray			*temp;
Daniel@0: 	
Daniel@0: 	zip = zi;
Daniel@0: 	ziend = zi + numElements;
Daniel@0: 	while (zip < ziend)
Daniel@0: 	{
Daniel@0: 		allImZero = allImZero && (*zip++ == 0.0);
Daniel@0: 		if (!allImZero)
Daniel@0: 			return;
Daniel@0: 	}
Daniel@0: 	
Daniel@0: 	temp = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
Daniel@0: 	z1p = zr;
Daniel@0: 	z2p = mxGetPr(temp);
Daniel@0: 	zrend = z1p + numElements;
Daniel@0: 	while (z1p < zrend)
Daniel@0: 		*z2p++ = *z1p++;
Daniel@0: 	mxDestroyArray(plhs[0]);
Daniel@0: 	plhs[0] = temp;
Daniel@0: 	return;
Daniel@0: }