#include <u.h>
#include <libc.h>
#include <bio.h>

#include "apl9.h"

Rune primfuncnames[] = L"+-×÷*⍟⌹○!?|⌈⌊⊥⊤⊣⊢=≠≤<>≥≡≢∨∧⍲⍱↑↓⊂⊃⊆⌷⍋⍒⍳⍸∊⍷∪∩~,⍪⍴⌽⊖⍉⍎⍕∇";

fnmonad monadfunctiondefs[] = {
	fnSame, /* + */
	fnNegate, /* - */
	fnSign, /* × */
	fnRecip, /* ÷ */
	fnExponent, /* * */
	fnNaturalLog, /* ⍟ */
	0, /* ⌹ */
	fnPiTimes, /* ○ */
	0, /* ! */
	0, /* ? */
	fnMagnitude, /* | */
	fnCeiling, /* ⌈ */
	fnFloor, /* ⌊ */
	0, /* ⊥ */
	0, /* ⊤ */
	fnSame, /* ⊣ */
	fnSame, /* ⊢ */
	0, /* = */
	0, /* ≠ */
	0, /* ≤ */
	0, /* < */
	0, /* > */
	0, /* ≥ */
	fnDepth, /* ≡ */
	fnTally, /* ≢ */
	0, /* ∨ */
	0, /* ∧ */
	0, /* ⍲ */
	0, /* ⍱ */
	fnMix, /* ↑ */
	fnSplit, /* ↓ */
	fnEnclose, /* ⊂ */
	fnDisclose, /* ⊃ */
	fnNest, /* ⊆ */
	0, /* ⌷ */
	fnGradeUp, /* ⍋ */
	fnGradeDown, /* ⍒ */
	fnIndexGenerator, /* ⍳ */
	0, /* ⍸ */
	fnEnlist, /* ∊ */
	0, /* ⍷ */
	0, /* ∪ */
	0, /* ∩ */
	fnNot, /* ~ */
	fnRavel, /* , */
	fnTable, /* ⍪ */
	fnShape, /* ⍴ */
	fnReverseLast, /* ⌽ */
	fnReverseFirst, /* ⊖ */
	fnTranspose, /* ⍉ */
	fnExecute, /* ⍎ */
	0, /* ⍕ */
	fnSelfReference1, /* ∇ */
};

fndyad dyadfunctiondefs[] = {
	fnPlus, /* + */
	fnMinus, /* - */
	fnTimes, /* × */
	fnDivide, /* ÷ */
	fnPower, /* * */
	fnLogarithm, /* ⍟ */
	0, /* ⌹ */
	0, /* ○ */
	0, /* ! */
	0, /* ? */
	fnResidue, /* | */
	fnMaximum, /* ⌈ */
	fnMinimum, /* ⌊ */
	0, /* ⊥ */
	0, /* ⊤ */
	fnLeft, /* ⊣ */
	fnRight, /* ⊢ */
	fnEqual, /* = */
	fnNotEqual, /* ≠ */
	fnLessEqual, /* ≤ */
	fnLess, /* < */
	fnGreater, /* > */
	fnGreaterEqual, /* ≥ */
	fnMatch, /* ≡ */
	0, /* ≢ */
	fnOr, /* ∨ */
	fnAnd, /* ∧ */
	fnNand, /* ⍲ */
	fnNor, /* ⍱ */
	fnTake, /* ↑ */
	0, /* ↓ */
	0, /* ⊂ */
	0, /* ⊃ */
	0, /* ⊆ */
	fnIndex, /* ⌷ */
	0, /* ⍋ */
	0, /* ⍒ */
	0, /* ⍳ */
	0, /* ⍸ */
	fnMembership, /* ∊ */
	0, /* ⍷ */
	fnUnion, /* ∪ */
	0, /* ∩ */
	fnExcluding, /* ~ */
	fnCatenateLast, /* , */
	fnCatenateFirst, /* ⍪ */
	fnReshape, /* ⍴ */
	0, /* ⌽ */
	0, /* ⊖ */
	0, /* ⍉ */
	0, /* ⍎ */
	0, /* ⍕ */
	fnSelfReference2, /* ∇ */
};

vlong gcd_int(vlong, vlong);
double gcd_float(double, double);

/* Runner function */
Array *
runfunc(Function f, Array *left, Array *right)
{
	Array *result;
	if(f.type == FunctypeDfn || (f.type == FunctypeOp && f.operator.type == OperatortypeDop)){
		Rune *code;
		Datum *lefto = nil;
		Datum *righto = nil;
		if(f.type == FunctypeDfn)
			code = f.dfn;
		else{
			lefto = f.operator.left;
			righto = f.operator.right;
			code = f.operator.dop;
		}

		pushdfnframe(code, lefto, righto, left, right);
		Datum *dfnres = evalline(code, nil, 0);
		popdfnframe();

		result = (*dfnres).array; /* TODO what if the evaluation failed */
	}else if(f.type == FunctypePrim){
		if(left){
			fndyad d = dyadfunctiondefs[f.code];
			if(d == nil){
				Rune *err = runesmprint("dyadic %C", primfuncnames[f.code]);
				throwerror(err, ENotImplemented);
			}
			result = d(left, right);
		}else{
			fnmonad m = monadfunctiondefs[f.code];
			if(m == nil){
				Rune *err = runesmprint("monadic %C", primfuncnames[f.code]);
				throwerror(err, ENotImplemented);
			}
			result = m(right);
		}
	}else if(f.type == FunctypeOp && f.operator.type == OperatortypeHybrid){
		opmonad m = hybridoperatordefs[f.operator.code];
		if(m == nil){
			Rune *err = runesmprint("monadic %C", primhybridnames[f.operator.code]);
			throwerror(err, ENotImplemented);
		}
		result = m(f.operator.left, left, right);
	}else if(f.type == FunctypeOp){
		/* TODO assumes prim op, not dop */
		if(f.operator.dyadic){
			opdyad d = dyadoperatordefs[f.operator.code];
			if(d == nil){
				Rune *err = runesmprint("dyadic %C", primdyadopnames[f.operator.code]);
				throwerror(err, ENotImplemented);
			}
			result = d(f.operator.left, f.operator.right, left, right);
		}else{
			opmonad m = monadoperatordefs[f.operator.code];
			if(m == nil){
				Rune *err = runesmprint("monadic %C", primmonopnames[f.operator.code]);
				throwerror(err, ENotImplemented);
			}
			result = m(f.operator.left, left, right);
		}
	}else if(f.type == FunctypeQuad){
		if(left){
			if(f.quad->dyadfn == nil){
				Rune *err = runesmprint("dyadic %S", f.quad->name);
				throwerror(err, ENotImplemented);
			}
			result = f.quad->dyadfn(left, right);
		}else{
			if(f.quad->monadfn == nil){
				Rune *err = runesmprint("monadic %S", f.quad->name);
				throwerror(err, ENotImplemented);
			}
			result = f.quad->monadfn(right);
		}
	}else if(f.type == FunctypeHybrid){
		fndyad d = hybridfunctiondefs[f.code];
		if(d == nil){
			Rune *err = runesmprint("dyadic %C", primhybridnames[f.code]);
			throwerror(err, ENotImplemented);
		}
		result = d(left, right);
	}else if(f.type == FunctypeTrain)
		result = runtrain(f.train.funcs, f.train.nfuncs, left, right, nil);
	else{
		throwerror(L"runfunc for this function type", ENotImplemented);
		result = nil;
	}
	Array *tmp = result;
	result = simplifyarray(result);
	freearray(tmp);
	return result;
}

Array *
rundfn(Rune *code, Datum *lefto, Datum *righto, Array *left, Array *right)
{
	if(lefto == nil && righto == nil){
		Function dfn;
		dfn.type = FunctypeDfn;
		dfn.dfn = code;
		return runfunc(dfn, left, right);
	}else if(lefto != nil){
		Function dop;
		dop.type = FunctypeOp;
		dop.operator.type = OperatortypeDop;
		dop.operator.left = lefto;
		dop.operator.right = righto;
		dop.operator.dyadic = righto != nil;
		dop.operator.dop = code;
		return runfunc(dop, left, right);
	}else{
		throwerror(L"Malformed call to rundfn", ENotImplemented);
		return nil;
	}
}

Array *
runtrain(Function *funcs, int nfuncs, Array *left, Array *right, Array *acc)
{
	if(nfuncs >= 3 && acc == nil && !(funcs[0].left != nil && nfuncs == 3)){
		Array *r = runfunc(funcs[nfuncs-1], left, right);
		Array *l = runfunc(funcs[nfuncs-3], left, right);
		Array *c = runfunc(funcs[nfuncs-2], l, r);
		freearray(r);
		freearray(l);
		if(nfuncs == 3)
			return c;
		else
			return runtrain(funcs, nfuncs-3, left, right, c);
	}else if(nfuncs >= 2 && acc != nil && !(funcs[0].left != nil && nfuncs == 2)){
		Array *l = runfunc(funcs[nfuncs-2], left, right);
		Array *c = runfunc(funcs[nfuncs-1], l, acc);
		freearray(acc);
		freearray(l);
		if(nfuncs == 2)
			return c;
		else
			return runtrain(funcs, nfuncs-2, left, right, c);
	}else if(nfuncs == 2 && acc == nil){
		Array *r = runfunc(funcs[1], left, right);
		Array *c = runfunc(funcs[0], funcs[0].left, r);
		freearray(r);
		return c;
	}else if(nfuncs == 1 && acc != nil){
		Array *c = runfunc(funcs[0], funcs[0].left, acc);
		freearray(acc);
		return c;
	}else{
		throwerror(L"train combination", ESyntax);
		return nil;
	}
}

/* Monadic functions */

Array *
fnNegate(Array *right)
{
	return rundfn(L"0-⍵", nil, nil, nil, right);
}

Array *
fnSign(Array *right)
{
	return rundfn(L"(¯1×⍵<0)+⍵>0", nil, nil, nil, right);
}

Array *
fnRecip(Array *right)
{
	return rundfn(L"1÷⍵", nil, nil, nil, right);
}

Array *
fnExponent(Array *right)
{
	Array *e = mkscalarfloat(2.7182818284590452353602);
	Array *res = fnPower(e, right);
	freearray(e);
	return res;
}

Array *
fnNaturalLog(Array *right)
{
	Array *e = mkscalarfloat(2.7182818284590452353602);
	Array *res = fnLogarithm(e, right);
	freearray(e);
	return res;
}

Array *
fnPiTimes(Array *right)
{
	Array *pi = mkscalarfloat(PI);
	Array *res = fnTimes(pi, right);
	freearray(pi);
	return res;
}

Array *
fnMagnitude(Array *right)
{
	return rundfn(L"⍵××⍵", nil, nil, nil, right);
}

Array *
fnCeiling(Array *right)
{
	return rundfn(L"-⌊-⍵", nil, nil, nil, right);
}

Array *
fnFloor(Array *right)
{
	Array *res = nil;
	switch(right->type){
	case AtypeInt:
		res = fnSame(right);
		break;
	case AtypeFloat:
		res = duparrayshape(right, AtypeInt);
		for(int i = 0; i < right->size; i++)
			res->intdata[i] = floor(right->floatdata[i]);
		break;
	case AtypeArray:
		res = duparrayshape(right, AtypeArray);
		for(int i = 0; i < right->size; i++)
			res->arraydata[i] = fnFloor(right->arraydata[i]);
		break;
	default:
		throwerror(nil, EType);
	}
	return res;
}

Array *
fnSame(Array *right)
{
	incref(right);
	return right;
}

Array *
fnDepth(Array *right)
{
	int uniform;
	int depth = arraydepth(right, &uniform);
	if(!uniform)
		depth = -depth;
	return mkscalarint(depth);
}

Array *
fnTally(Array *right)
{
	return mkscalarint(right->rank==0 ? 1 : right->shape[0]);
}

Array *
fnMix(Array *right)
{
	if(right->type != AtypeArray || right->size == 0)
		return fnSame(right);

	int commonrank = 0;
	int i,j;
	for(i = 0; i < right->size; i++)
		if(right->arraydata[i]->rank > commonrank)
			commonrank = right->arraydata[i]->rank;

	Array *commonshape = allocarray(AtypeInt, 1, commonrank);
	commonshape->shape[0] = commonrank;
	for(i = 0; i < commonrank; i++)
		commonshape->intdata[i] = 0;

	for(i = 0; i < right->size; i++){
		Array *a = right->arraydata[i];
		for(j = 0; j < a->rank; j++){
			if(a->shape[a->rank-1-j] > commonshape->intdata[commonrank-1-j])
				commonshape->intdata[commonrank-1-j] = a->shape[a->rank-1-j];
		}
	}

	int size = 1;
	int commonsize = 1;
	for(i = 0; i < right->rank; i++)
		size *= right->shape[i];
	for(i = 0; i < commonshape->size; i++){
		size *= commonshape->intdata[i];
		commonsize *= commonshape->intdata[i];
	}

	/* TODO: think about types */
	Array *result = allocarray(right->arraydata[0]->type, right->rank + commonrank, size);
	for(i = 0; i < right->rank; i++)
		result->shape[i] = right->shape[i];
	for(j = 0; j < commonshape->size; j++)
		result->shape[i+j] = commonshape->intdata[j];

	int *index = malloc(sizeof(int) * commonrank);
	int offset;
	for(i = 0; i < size/commonsize; i++){
		Array *a = right->arraydata[i];
		Array *fill = fillelement(a);
		if(a->rank == 0){
			memcpy(
				result->rawdata + i * commonsize * datasizes[a->type],
				a->rawdata, datasizes[a->type]);
			if(a->type == AtypeArray)
				incref(a->arraydata[0]);
			for(j = 1; j < commonsize; j++){
				memcpy(result->rawdata + (i * commonsize + j) * datasizes[a->type],
					fill->rawdata, datasizes[a->type]);
				if(fill->type == AtypeArray)
					incref(fill->arraydata[0]);
			}
		}else{
			for(j = 0; j < commonrank; j++)
				index[j] = 0;
			for(j = 0, offset = 0; offset < commonsize; j++){
				for(int k = 0; index[commonrank-1-k] == a->shape[a->rank-1-k]; k++){
					int nfill = commonshape->intdata[commonrank-1-k] - a->shape[a->rank-1-k];
					if(nfill)
						print("Adding %d fills\n", nfill);
					while(nfill--){
						memcpy(result->rawdata + (i * commonsize + offset) * datasizes[a->type],
							fill->rawdata, datasizes[a->type]);
						if(fill->type == AtypeArray)
							incref(fill->arraydata[0]);
						offset++;
					}
					index[commonrank-1-k] = 0;
					index[commonrank-2-k]++;
				}
				if(offset < commonsize){
					memcpy(
						result->rawdata + (i * commonsize + offset) * datasizes[a->type],
						a->rawdata + j * datasizes[a->type], datasizes[a->type]);
					if(a->type == AtypeArray)
						incref(a->arraydata[j]);
					offset++;
					index[commonrank-1]++;
				}
			}
		}
		freearray(fill);
	}
	free(index);
	free(commonshape);
	return result;
}

Array *
fnSplit(Array *right)
{
	Rune *code = L"0≡≢⍴⍵: ⍵	⋄ 1≡≢⍴⍵: ⊂⍵ ⋄ (⊂⍵)⌷⍨¨⍳≢⍵";	
	return rundfn(code, nil, nil, nil, right);
}

Array *
fnEnclose(Array *right)
{
	incref(right);
	if(simplescalar(right))
		return right;
	else{
		Array *res = allocarray(AtypeArray, 0, 1);
		res->arraydata[0] = right;
		return res;
	}
}

Array *
fnDisclose(Array *right)
{
	if(right->size == 0)
		return fillelement(right);
	else
		return arrayitem(right, 0);
}

Array *
fnNest(Array *right)
{
	if(simplearray(right))
		return fnEnclose(right);
	else
		return fnSame(right);
}

Array *
fnGradeUp(Array *right)
{
	if(right->rank == 0)
		throwerror(nil, ERank);

	int i,j;
	int len = right->shape[0];
	Array **elems = malloc(sizeof(Array *) * len);
	Array *index = mkscalarint(globalIO());
	Array *order = allocarray(AtypeInt, 1, len);
	order->shape[0] = len;

	vlong io = globalIO();
	for(i = 0; i < len; i++, index->intdata[0]++){
		order->intdata[i] = io + i;
		elems[i] = fnIndex(index, right);
	}

	/* Do a insertion sort on elems, while also updating order */
	for(i = 1; i < len; i++){
		for(j = i; j > 0 && comparearray(elems[j], elems[j-1], 0) == -1; j--){
			Array *tmpA = elems[j];
			elems[j] = elems[j-1];
			elems[j-1] = tmpA;

			int tmpI = order->intdata[j];
			order->intdata[j] = order->intdata[j-1];
			order->intdata[j-1] = tmpI;
		}
	}

	freearray(index);
	for(i = 0; i < len; i++)
		freearray(elems[i]);
	free(elems);
	return order;
}

Array *
fnGradeDown(Array *right)
{
	Array *tmp = fnGradeUp(right);
	Array *res = fnReverseFirst(tmp);
	freearray(tmp);
	return res;
}

Array *
fnIndexGenerator(Array *right)
{
	/* TODO only works for creating vectors */
	vlong n = right->intdata[0];
	Array *res = allocarray(AtypeInt, 1, n);
	res->shape[0] = n;
	vlong io = globalIO();
	for(vlong i = 0; i < n; i++)
		res->intdata[i] = i + io;
	return res;
}

Array *
fnEnlist(Array *right)
{
	if(right->size == 0)
		return fnSame(right);
	if(right->type != AtypeArray)
		return fnRavel(right);
	else{
		Array *res = fnEnlist(right->arraydata[0]);
		for(int i = 1; i < right->size; i++){
			Array *old = res;
			Array *tmp = fnEnlist(right->arraydata[i]);
			res = fnCatenateFirst(res, tmp);
			freearray(old);
			freearray(tmp);
		}
		return res;
	}
}

Array *
fnNot(Array *right)
{
	if(right->type != AtypeInt)
		throwerror(nil, EType);
	Array *res = duparray(right);
	for(int i = 0; i < res->size; i++){
		if(res->intdata[i] == 0)
			res->intdata[i] = 1;
		else if(res->intdata[i] == 1)
			res->intdata[i] = 0;
		else{
			freearray(res);
			throwerror(nil, EType);
		}
	}
	return res;
}

Array *
fnRavel(Array *right)
{
	Array *res = duparray(right);
	res->rank = 1;
	res->shape = realloc(res->shape, sizeof(int) * 1);
	res->shape[0] = res->size;
	return res;
}

Array *
fnTable(Array *right)
{
	Array *res = duparray(right);
	res->rank = 2;
	res->shape = realloc(res->shape, sizeof(int) * 2);
	res->shape[0] = right->rank ? res->shape[0] : 1;
	res->shape[1] = right->size / res->shape[0];
	return res;
}

Array *
fnShape(Array *right)
{
	Array *res = allocarray(AtypeInt, 1, right->rank);
	res->shape[0] = right->rank;
	for(int i = 0; i < right->rank; i++)
		res->intdata[i] = right->shape[i];
	return res;
}

Array *
fnReverseLast(Array *right)
{
	if(right->rank < 1)
		return fnSame(right);

	Array *res = duparray(right);

	int nrows = 1;
	int rowsize = res->shape[res->rank-1];
	for(int i = 0; i < res->rank - 1; i++)
		nrows *= res->shape[i];

	for(int row = 0; row < nrows; row++){
		for(int i = 0; i < rowsize; i++)
			memcpy(
				res->rawdata + (row * rowsize + i) * datasizes[res->type],
				right->rawdata + ((1+row) * rowsize - 1 - i) * datasizes[res->type],
				datasizes[res->type]);
	}
	return res;
}

Array *
fnReverseFirst(Array *right)
{
	if(right->rank < 1 || right->shape[0] == 0)
		return fnSame(right);

	Array *res = duparray(right);
	int cells = res->shape[0];
	int elems = res->size / cells;
	for(int i = 0; i < cells; i++)
		memcpy(
			res->rawdata + i * elems * datasizes[res->type],
			right->rawdata + (cells - 1 - i) * elems * datasizes[res->type],
			datasizes[res->type] * elems);
	return res;
}

Array *
fnTranspose(Array *right)
{
	Array *res = duparray(right);
	for(int i = 0; i < res->rank; i++)
		res->shape[i] = right->shape[res->rank - 1 - i];

	int from, to;
	int *sizesFrom = malloc(sizeof(int) * right->rank);
	int *sizesTo = malloc(sizeof(int) * right->rank);
	int accFrom = 1, accTo = 1;
	for(int i = 0; i < right->rank; i++){
		sizesFrom[i] = accFrom;
		sizesTo[i] = accTo;
		accFrom *= res->shape[i];
		accTo *= right->shape[i];
	}

	for(from = 0; from < right->size; from++){
		to = 0;
		int tmp = from;
		for(int i = right->rank-1; i >= 0; i--){
			to += sizesTo[right->rank-1-i] * (tmp / sizesFrom[i]);
			tmp = tmp % sizesFrom[i];
		}
		memcpy(
			res->rawdata + to * datasizes[res->type],
			right->rawdata + from * datasizes[res->type],
			datasizes[res->type]);
	}
	free(sizesFrom);
	free(sizesTo);
	return res;
}

Array *
fnExecute(Array *right)
{
	if(right->type != AtypeRune)
		throwerror(nil, EType);
	Rune *code = pparray(right);
	Datum *result = evalline(code, nil, 1);
	free(code);
	if(!result)
		throwerror(L"No result produced by ⍎", EDomain);
	if(result->tag != ArrayTag)
		throwerror(L"Result of ⍎ must be an array", EType);
	Array *res = result->array;
	free(result);
	return res;
}

Array *
fnSelfReference1(Array *right)
{
	DfnFrame *dfn = getcurrentdfn();
	if(dfn)
		return rundfn(dfn->code, dfn->lefto, dfn->righto, nil, right);
	else{
		throwerror(nil, ESyntax);
		return nil;
	}
}

/* Dyadic functions */

/* macro to define dyadic scalar functions */
#define SCALAR_FUNCTION_2(name, forcefloat, restype, cases) \
Array *name(Array *left, Array *right){\
	Array *leftarr, *rightarr;\
	if(!commontype(left, right, &leftarr, &rightarr, forcefloat)) throwerror(nil, EType);\
	if(!scalarextend(leftarr, rightarr, &left, &right)) throwerror(L"scalar extension fail", ERank);\
	Array *res;\
	if(left->type != AtypeArray && restype != left->type)\
		res = duparrayshape(left, restype);\
	else\
		res = duparray(left);\
	for(int i = 0; i < left->size; i++)\
		switch(left->type){\
		default: throwerror(nil, EType); break;\
		case AtypeArray:\
			freearray(res->arraydata[i]);\
			res->arraydata[i] = name(left->arraydata[i], right->arraydata[i]);\
			break;\
		cases\
		}\
	freearray(leftarr); freearray(rightarr); freearray(left); freearray(right);\
	return res;}

SCALAR_FUNCTION_2(fnPlus, 0, left->type,
	case AtypeFloat:
		res->floatdata[i] += right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] += right->intdata[i];
		break;
)

SCALAR_FUNCTION_2(fnMinus, 0, left->type,
	case AtypeFloat:
		res->floatdata[i] -= right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] -= right->intdata[i];
		break;
)

SCALAR_FUNCTION_2(fnTimes, 0, left->type,
	case AtypeFloat:
		res->floatdata[i] *= right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] *= right->intdata[i];
		break;
)

SCALAR_FUNCTION_2(fnDivide, 1, left->type,
	case AtypeFloat:
		if(right->floatdata[i] == 0){
			if(globalDIV())
				res->floatdata[i] = 0;
			else
				throwerror(nil, EDomain);
		}
		else
			res->floatdata[i] /= right->floatdata[i];
		break;
)

SCALAR_FUNCTION_2(fnPower, 0, left->type,
	case AtypeFloat:
		res->floatdata[i] = pow(res->floatdata[i], right->floatdata[i]);
		break;
	case AtypeInt:
		res->intdata[i] = pow(res->intdata[i], right->intdata[i]);
		break;
)

SCALAR_FUNCTION_2(fnLogarithm, 1, left->type,
	case AtypeFloat:
		res->floatdata[i] = log(right->floatdata[i])/log(res->floatdata[i]);
		break;
)

SCALAR_FUNCTION_2(fnResidue, 1, left->type,
	case AtypeFloat:
		if(res->floatdata[i] == 0)
			res->floatdata[i] = right->floatdata[i];
		else
			res->floatdata[i] = right->floatdata[i] - res->floatdata[i] * floor(right->floatdata[i]/res->floatdata[i]);
)

SCALAR_FUNCTION_2(fnMaximum, 0, left->type,
	case AtypeFloat:
		if(res->floatdata[i] < right->floatdata[i])
			res->floatdata[i] = right->floatdata[i];
	case AtypeInt:
		if(res->intdata[i] < right->intdata[i])
			res->intdata[i] = right->intdata[i];		
)

SCALAR_FUNCTION_2(fnMinimum, 0, left->type,
	case AtypeFloat:
		if(res->floatdata[i] > right->floatdata[i])
			res->floatdata[i] = right->floatdata[i];
	case AtypeInt:
		if(res->intdata[i] > right->intdata[i])
			res->intdata[i] = right->intdata[i];		
)

Array *
fnLeft(Array *left, Array *right)
{
	USED(right);
	incref(left);
	return left;
}

Array *
fnRight(Array *left, Array *right)
{
	USED(left);
	incref(right);
	return right;
}

SCALAR_FUNCTION_2(fnEqual, 0, AtypeInt,
	case AtypeFloat:
		res->intdata[i] = left->floatdata[i] == right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] = left->intdata[i] == right->intdata[i];
		break;
	case AtypeRune:
		res->intdata[i] = left->runedata[i] == right->runedata[i];
		break;
)

SCALAR_FUNCTION_2(fnNotEqual, 0, AtypeInt,
	case AtypeFloat:
		res->intdata[i] = left->floatdata[i] != right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] = left->intdata[i] != right->intdata[i];
		break;
	case AtypeRune:
		res->intdata[i] = left->runedata[i] != right->runedata[i];
		break;
)

SCALAR_FUNCTION_2(fnLessEqual, 0, AtypeInt,
	case AtypeFloat:
		res->intdata[i] = left->floatdata[i] <= right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] = left->intdata[i] <= right->intdata[i];
		break;
)

SCALAR_FUNCTION_2(fnLess, 0, AtypeInt,
	case AtypeFloat:
		res->intdata[i] = left->floatdata[i] < right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] = left->intdata[i] < right->intdata[i];
		break;
)

SCALAR_FUNCTION_2(fnGreater, 0, AtypeInt,
	case AtypeFloat:
		res->intdata[i] = left->floatdata[i] > right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] = left->intdata[i] > right->intdata[i];
		break;
)

SCALAR_FUNCTION_2(fnGreaterEqual, 0, AtypeInt,
	case AtypeFloat:
		res->intdata[i] = left->floatdata[i] >= right->floatdata[i];
		break;
	case AtypeInt:
		res->intdata[i] = left->intdata[i] >= right->intdata[i];
		break;
)

Array *
fnMatch(Array *left, Array *right)
{
	int cmp = comparearray(left, right, 1);
	return mkscalarint(cmp == 0);
}

SCALAR_FUNCTION_2(fnOr, 0, left->type,
	case AtypeInt:
		if((left->intdata[i] == 0 || left->intdata[i] == 1) && (right->intdata[i] == 0 || right->intdata[i] == 1))
			res->intdata[i] = left->intdata[i] || right->intdata[i];
		else
			res->intdata[i] = gcd_int(left->intdata[i], right->intdata[i]);
		break;
	case AtypeFloat:
		res->floatdata[i] = gcd_float(left->floatdata[i], right->floatdata[i]);
		break;
)

SCALAR_FUNCTION_2(fnAnd, 0, left->type,
	case AtypeInt:
		if((left->intdata[i] == 0 || left->intdata[i] == 1) && (right->intdata[i] == 0 || right->intdata[i] == 1))
			res->intdata[i] = left->intdata[i] && right->intdata[i];
		else{
			vlong p = left->intdata[i] * right->intdata[i];
			if(p < 0)
				p = -p;
			res->intdata[i] = p / gcd_int(left->intdata[i], right->intdata[i]);
		}
		break;
	case AtypeFloat:
		res->floatdata[i] = fabs(left->floatdata[i] * right->floatdata[i]) / gcd_float(left->floatdata[i], right->floatdata[i]);
		break;
)

SCALAR_FUNCTION_2(fnNand, 0, AtypeInt,
	case AtypeInt:
		if(left->intdata[i] != 0 && left->intdata[i] != 1)
			throwerror(nil, EType);
		if(right->intdata[i] != 0 && right->intdata[i] != 1)
			throwerror(nil, EType);
		res->intdata[i] = !(left->intdata[i] && right->intdata[i]);
		break;
)

SCALAR_FUNCTION_2(fnNor, 0, AtypeInt,
	case AtypeInt:
		if(left->intdata[i] != 0 && left->intdata[i] != 1)
			throwerror(nil, EType);
		if(right->intdata[i] != 0 && right->intdata[i] != 1)
			throwerror(nil, EType);
		res->intdata[i] = !(left->intdata[i] || right->intdata[i]);
		break;
)

Array *
fnTake(Array *left, Array *right)
{
	if(left->type != AtypeInt)
		throwerror(nil, EType);
	if(left->rank > 1)
		throwerror(nil, ERank);
	if(left->size > right->rank)
		throwerror(nil, ELength);

	int i;
	if(left->size == right->rank)
		left = fnSame(left);
	else{
		Array *old = left;
		left = fnShape(right);
		for(i = 0; i < old->size; i++)
			left->intdata[i] = old->intdata[i];
	}

	if(right->rank == 0){
		Array *leftshape = fnShape(left);
		right = fnReshape(leftshape, right);
		freearray(leftshape);
	}else
		right = fnSame(right);
	
	int *shape = malloc(sizeof(int) * left->size);
	int size = 1;
	for(i = 0; i < left->size; i++){
		int s = left->intdata[i];
		shape[i] = s < 0 ? -s : s;
		size *= shape[i];
	}

	Array *result = allocarray(right->type, right->rank, size);
	for(i = 0; i < right->rank; i++)
		result->shape[i] = shape[i];

	int *index = mallocz(sizeof(int) * left->size, 1);
	int offset;
	for(i = 0, offset = 0; offset < size; i++){
		for(int j = left->size-1; index[j] == shape[j]; j--){
			index[j] = 0;
			index[j-1]++;
		}
		print("Result Index: ");
		for(int j = 0; j < left->size; j++)
			print("%d ", index[j]);
		print("\n");

		/* if index is part of left vector, select those places */

		offset++;
		index[left->size-1]++;
	}

	freearray(left);
	freearray(right);
	return result;
}

Array *
fnIndex(Array *left, Array *right)
{
	int io = globalIO();
	int i;

	if(left->rank > 1)
		throwerror(nil, ERank);
	if(left->type != AtypeArray && left->type != AtypeInt)
		throwerror(nil, EType);
	if(left->size > right->rank)
		throwerror(nil, ELength);

	/* extend left index vector to full format */
	Array *oldleft = left;
	left = allocarray(AtypeArray, 1, right->rank);
	left->shape[0] = right->rank;
	for(i = 0; i < left->size; i++){
		if(i >= oldleft->size){
			Array *n = mkscalarint(right->shape[i]);
			left->arraydata[i] = fnIndexGenerator(n);
			freearray(n);
		}else if(oldleft->type == AtypeInt){
			if(oldleft->intdata[i] < io || oldleft->intdata[i] >= io + right->shape[i])
				throwerror(nil, EIndex);
			left->arraydata[i] = mkscalarint(oldleft->intdata[i]);
		}else if(oldleft->type == AtypeArray){
			Array *sub = oldleft->arraydata[i];
			if(sub->type != AtypeInt)
				throwerror(nil, EType);
			for(int j = 0; j < sub->size; j++){
				if(sub->intdata[j] < io || sub->intdata[j] >= io + right->shape[i])
					throwerror(nil, EIndex);
			}
			left->arraydata[i] = oldleft->arraydata[i];
			incref(left->arraydata[i]);
		}
	}

	Array *shape = allocarray(AtypeInt, 1, 0);
	shape->shape[0] = 0;

	int rank = 0;
	int size = 1;
	for(i = 0; i < left->size; i++){
		if(left->type == AtypeArray){
			Array *tmp = shape;
			Array *newShape = fnShape(left->arraydata[i]);
			shape = fnCatenateFirst(tmp, newShape);
			freearray(tmp);
			freearray(newShape);
		}
	}

	for(i = 0; i < shape->size; i++){
		size *= shape->intdata[i];
		rank++;
	}

	Array *result = allocarray(right->type, rank, size);
	for(i = 0; i < rank; i++)
		result->shape[i] = shape->intdata[i];
	freearray(shape);
	int *leftindex = mallocz(sizeof(int) * right->rank, 1);
	for(i = 0; i < result->size; i++){
		for(int j = left->size-1; leftindex[j] == left->arraydata[j]->size; j--){
			leftindex[j] = 0;
			leftindex[j-1]++;
		}
		int from = 0;
		int cellsize = right->size;
		for(int j = 0; j < left->size; j++){
			cellsize = cellsize / right->shape[j];
			from += cellsize * (left->arraydata[j]->intdata[leftindex[j]] - io);
		}
		memcpy(
			result->rawdata + i * datasizes[result->type],
			right->rawdata + from * datasizes[result->type],
			datasizes[result->type]);
		if(result->type == AtypeArray)
			incref(result->arraydata[i]);
		leftindex[left->size-1]++;
	}
	free(leftindex);
	freearray(left);
	return result;
}

Array *
fnMembership(Array *left, Array *right)
{
	/* TODO avoid outer product */
	return rundfn(L"∨/⍺≡⌾⍵", nil, nil, left, right);
}

Array *
fnUnion(Array *left, Array *right)
{
	if(left->rank > 1 || right->rank > 1)
		throwerror(nil, ERank);
	return rundfn(L"⍺⍪⍵~⍺", nil, nil, left, right);
}

Array *
fnExcluding(Array *left, Array *right)
{
	return rundfn(L"(~⍺∊⍵)⌿⍺", nil, nil, left, right);
}

Array *
fnCatenateLast(Array *left, Array *right)
{
	return rundfn(L"⍉(⍉⍺)⍪⍉⍵", nil, nil, left, right);
}

Array *
fnCatenateFirst(Array *left, Array *right)
{
	Array *leftarr = nil;
	Array *rightarr = nil;

	if(left->rank == 0 && right->rank != 0){
		/* extend left to right->rank with first axis=1 */
		rightarr = fnSame(right);
		Array *shape = fnShape(right);
		shape->intdata[0] = 1;
		leftarr = fnReshape(shape, left);
		freearray(shape);
	}else if(left->rank != 0 && right->rank == 0){
		/* extend right to left->rank with first axis=1 */
		leftarr = fnSame(left);
		Array *shape = fnShape(left);
		shape->intdata[0] = 1;
		rightarr = fnReshape(shape, right);
		freearray(shape);
	}else if(left->rank == 0 && right->rank == 0){
		/* turn both scalars into vectors */
		leftarr = fnRavel(left);
		rightarr = fnRavel(right);
	}else{
		/* Check that the shapes match */
		if(left->rank == right->rank-1){
			/* extend left with unit dimension */
			Array *shape = allocarray(AtypeInt, 1, left->rank+1);
			shape->intdata[0] = 1;
			for(int i = 1; i < left->rank+1; i++)
				shape->intdata[i] = left->shape[i-1];
			leftarr = fnReshape(shape, left);
			rightarr = fnSame(right);
			freearray(shape);
		}else if(right->rank == left->rank-1){
			/* extend right with unit dimension */
			Array *shape = allocarray(AtypeInt, 1, right->rank+1);
			shape->intdata[0] = 1;
			for(int i = 1; i < right->rank+1; i++)
				shape->intdata[i] = right->shape[i-1];
			rightarr = fnReshape(shape, right);
			leftarr = fnSame(left);
			freearray(shape);
		}else if(right->rank == left->rank){
			leftarr = fnSame(left);
			rightarr = fnSame(right);
		}else
			throwerror(nil, ERank);

		for(int i = 1; i < leftarr->rank; i++)
			if(leftarr->shape[i] != rightarr->shape[i])
				throwerror(nil, EShape);
	}

	int type, rank, leftsize, rightsize;
	if(leftarr->type == AtypeArray || rightarr->type == AtypeArray || leftarr->type != rightarr->type)
		type = AtypeArray;
	else
		type = leftarr->type;
	if(leftarr->rank > rightarr->rank)
		rank = leftarr->rank;
	else
		rank = rightarr->rank;
	
	leftsize = leftarr->shape[0];
	rightsize = rightarr->shape[0];
	for(int i = 1; i < rank; i++)
		leftsize *= leftarr->shape[i];
	for(int i = 1; i < rank; i++)
		rightsize *= rightarr->shape[i];

	Array *result = allocarray(type, rank, leftsize + rightsize);
	int i, j;
	result->shape[0] = leftarr->shape[0] + rightarr->shape[0];
	for(i = 1; i < result->rank; i++)
		result->shape[i] = leftarr->shape[i];

	/* TODO reduce duplicated code between copies from left and right */
	/* Copy data from the left array */
	for(i = 0, j = 0; i < leftarr->size; i++, j++){
		if(type == AtypeArray && leftarr->type == AtypeArray){
			result->arraydata[j] = leftarr->arraydata[i];
			incref(result->arraydata[j]);
		}else if(type == AtypeArray && leftarr->type != AtypeArray){
			result->arraydata[j] = arrayitem(leftarr, i);
		}else{
			memcpy(
				result->rawdata + j * datasizes[type],
				leftarr->rawdata + i * datasizes[type],
				datasizes[type]);
		}
	}

	/* Copy data from the right array */
	for(i = 0; i < rightarr->size; i++, j++){
		if(type == AtypeArray && rightarr->type == AtypeArray){
			result->arraydata[j] = rightarr->arraydata[i];
			incref(result->arraydata[j]);
		}else if(type == AtypeArray && rightarr->type != AtypeArray){
			result->arraydata[j] = arrayitem(rightarr, i);
		}else{
			memcpy(
				result->rawdata + j * datasizes[type],
				rightarr->rawdata + i * datasizes[type],
				datasizes[type]);
		}
	}

	freearray(leftarr);
	freearray(rightarr);
	return result;
}

Array *
fnReshape(Array *left, Array *right)
{
	vlong size = 1;
	int i;
	char *p;
	for(i = 0; i < left->size; i++)
		size *= left->intdata[i];
	if(left->size == 0)
		return arrayitem(right, 0);

	Array *res = allocarray(right->type, left->size, size);
	for(i = 0; i < left->size; i++)
		res->shape[i] = left->intdata[i];
	for(i = 0, p = res->rawdata; i < size; i++, p += datasizes[res->type])
		memcpy(p, right->rawdata + (datasizes[res->type] * (i % right->size)), datasizes[res->type]);
	if(res->type == AtypeArray)
		for(i = 0; i < res->size; i++)
			incref(res->arraydata[i]);
	return res;
}

Array *
fnSelfReference2(Array *left, Array *right)
{
	DfnFrame *dfn = getcurrentdfn();
	if(dfn)
		return rundfn(dfn->code, dfn->lefto, dfn->righto, left, right);
	else{
		throwerror(nil, ESyntax);
		return nil;
	}
}

/* helper functions */
vlong
gcd_int(vlong a, vlong b)
{
	if(b == 0)
		return a;
	else
		return gcd_int(b, a % b);
}

double
gcd_float(double a, double b)
{
	if(b < 0.0000000001 && b > -0.0000000001) /* Some error tolerance */
		return a;
	else
		return gcd_float(b, fmod(a,b));
}