/* rep_mult.c  repmat first two operands to the size provided by */
/* the third operand, then perform point multiply                */
/* 3 input, 1 output                                             */
/* C = rep_mult(A, B, sizes)                                     */

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	double		*xp, *yp, *zp, *pSizes;
	int			xnd, ynd, numElements = 1;
	const int	*xdim, *ydim;
	int         i, j, ndim;
	int			*s, *sx, *sy, *cpsx, *cpsy;
	int			*subs, *s1, *cpsx2, *cpsy2;

	if (nrhs != 3)
		mexErrMsgTxt("Incorrect number of inputs.");
	
	if (nlhs > 1)
		mexErrMsgTxt("Too many output arguments.");
	
	xnd = mxGetNumberOfDimensions(prhs[0]);
	ynd = mxGetNumberOfDimensions(prhs[1]);
	xdim = mxGetDimensions(prhs[0]);
	ydim = mxGetDimensions(prhs[1]);
	ndim = mxGetNumberOfElements(prhs[2]);

	pSizes = mxGetPr(prhs[2]);

	sx = (int *)malloc(sizeof(int)*ndim);
	sy = (int *)malloc(sizeof(int)*ndim);
	s =  (int *)malloc(sizeof(int)*ndim);
	s1 = (int *)malloc(sizeof(int)*ndim);
	*(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1;
	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
	subs =   (int *)malloc(sizeof(int)*ndim);
	cpsx2 =  (int *)malloc(sizeof(int)*ndim);
	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
	for(i=0; i<ndim; i++){
		subs[i] = 0;
		sx[i] = (i < xnd) ? xdim[i] : 1;
		sy[i] = (i < ynd) ? ydim[i] : 1;
		s[i] = (int)pSizes[i];
		s1[i] = s[i] - 1;
		numElements *= s[i];
	}
				
	for(i=0; i<ndim-1; i++){
		cpsx[i+1] = cpsx[i]*sx[i]--;
		cpsy[i+1] = cpsy[i]*sy[i]--;
		cpsx2[i] = cpsx[i]*sx[i];
		cpsy2[i] = cpsy[i]*sy[i];
	}
	cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
	
	plhs[0] = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
	zp = mxGetPr(plhs[0]);
	xp = mxGetPr(prhs[0]);
	yp = mxGetPr(prhs[1]);

	for(j=0; j<numElements; j++){
		*zp++ = *xp * *yp;
		for(i=0; i<ndim; i++){
			if(subs[i] == s1[i]){
				subs[i] = 0;
				if(sx[i])
					xp -= cpsx2[i];
				if(sy[i])
					yp -= cpsy2[i];
			}
			else{
				subs[i]++;
				if(sx[i])
					xp += cpsx[i];
				if(sy[i])
					yp += cpsy[i];
				break;
			}
		}
	}
	free(sx);
	free(sy);
	free(s);
	free(s1);
	free(cpsx);
	free(cpsy);
	free(subs);
	free(cpsx2);
	free(cpsy2);
}