/*
 * Copyright 1995,96 Thierry Bousch
 * Licensed under the Gnu Public License, Version 2
 *
 * $Id: Tensor.c,v 2.1 1996/08/18 09:28:55 bousch Exp $
 *
 * Operations on tensors
 */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "saml.h"
#include "saml-errno.h"
#include "mnode.h"
#include "builtin.h"

typedef struct {
	s_mnode *lit;
	int start, len;
} index_range;

typedef struct {
	struct mnode_header hdr;
	int indices;
	index_range idr[0];
} tensor_header;

#define T_HEADER(n)	((tensor_header*)(n))
#define INDICES(n)	(T_HEADER(n)->indices)
#define IRANGE(n)	(T_HEADER(n)->idr)
#define T_DATA(n)	((s_mnode**)(IRANGE(n)+INDICES(n)))
#define SCALAR(n)	(*T_DATA(n))

static void tensor_free (s_mnode*);
static s_mnode* tensor_build (const char*);
static gr_string* tensor_stringify (s_mnode*);
static s_mnode* tensor_make (s_mnode*);
static s_mnode* tensor_add (s_mnode*, s_mnode*);
static s_mnode* tensor_mul (s_mnode*, s_mnode*);
static s_mnode* tensor_div (s_mnode*, s_mnode*);
static s_mnode* tensor_gcd (s_mnode*, s_mnode*);
static int tensor_notzero (s_mnode*);
static s_mnode* tensor_zero (s_mnode*);
static s_mnode* tensor_negate (s_mnode*);
static s_mnode* tensor_one (s_mnode*);
static s_mnode* tensor2tensor (s_mnode *mn1, s_mnode *model);

s_mnode* tensor_move_literal (s_mnode*, s_mnode*, s_mnode*);
s_mnode* tensor_subs (s_mnode*, s_mnode*, s_mnode*);
s_mnode* tensor_diff (s_mnode*, s_mnode*);
s_mnode* tensor_sylvester (s_mnode*, s_mnode*, s_mnode*);

static s_mtype MathType_Tensor = {
	"Tensor",
	tensor_free, tensor_build, tensor_stringify,
	tensor_make, NULL,
	tensor_add, mn_std_sub, tensor_mul, tensor_div, tensor_gcd,
	tensor_notzero, NULL, NULL, mn_std_differ, NULL,
	tensor_zero, tensor_negate, tensor_one, NULL, NULL
};

void init_MathType_Tensor (void)
{
	register_mtype(ST_TENSOR, &MathType_Tensor);
	register_CV_routine(ST_TENSOR, ST_TENSOR, tensor2tensor);
}

static void copy_index_range (s_mnode* mn, s_mnode* mn1)
{
	int nbind = INDICES(mn);
	index_range *range = IRANGE(mn), *range1 = IRANGE(mn1);

	if (mn != mn1)
		memcpy(range, range1, nbind * sizeof(index_range));
	while (nbind--)
		copy_mnode((range++)->lit);
}

static s_mnode** get_tensor_entry (s_mnode* mn, const int *array)
{
	int i, nbind = INDICES(mn), position = 0, offset;
	index_range *range = IRANGE(mn);

	for (i = 0; i < nbind; i++) {
		offset = (*array) - (range->start);
		if (offset < 0 || offset >= range->len)
			return NULL;
		position = position * (range->len) + offset;
		++range; ++array;
	}
	return T_DATA(mn) + position;
}

static s_mnode* create_tensor (int nbind, int nb_elts)
{
	size_t bytes = sizeof(tensor_header) + nbind * sizeof(index_range)
		+ nb_elts * sizeof(s_mnode*);
	s_mnode* mn = __mnalloc(ST_TENSOR, bytes);
	INDICES(mn) = nbind;
	return mn;
}

static int count_tensor_elements (s_mnode *tensor)
{
	int nb_elts = 1, nbind = INDICES(tensor);
	index_range *range = IRANGE(tensor);

	while (nbind--)
		nb_elts *= (*range++).len;
	return nb_elts;
}

static s_mnode* tensor_mul_scalar (s_mnode* mn1, s_mnode* lambda, int flag)
{
	int nb_elts, nbind = INDICES(mn1);
	s_mnode *mn, **data, **data1;

	if (!mnode_notzero(lambda)) {
	    /* Special case: multiplication or division by zero */
	    if (flag)
		return mnode_error(SE_DIVZERO, "tensor_mul_scalar");
	    else
	    	return tensor_zero(mn1);
	}
	nb_elts = count_tensor_elements(mn1);
	mn = create_tensor(nbind, nb_elts);
	copy_index_range(mn, mn1);
	data1 = T_DATA(mn1); data = T_DATA(mn);
	if (flag) {
		/* Division */
		while (nb_elts--)
			*data++ = mnode_div(*data1++, lambda);
	} else {
		/* Multiplication */
		while (nb_elts--)
			*data++ = mnode_mul(*data1++, lambda);
	}
	return mn;
}

static void tensor_free (s_mnode *mn)
{
	int i, nb_elts = 1, nbind = INDICES(mn);
	index_range *irlist = IRANGE(mn);
	s_mnode **data = T_DATA(mn);

	for (i = 0; i < nbind; i++,irlist++) {
		nb_elts *= irlist->len;
		unlink_mnode(irlist->lit);
	}
	while (nb_elts--)
		unlink_mnode(*data++);
	free(mn);
}

static s_mnode* tensor_build (const char *str)
{
	char *bang, *litname;
	s_mnode *mn, *s;
	int index, litlen;
	index_range *range;

	if ((bang = strrchr(str, '!')) == NULL) {
		/* It is a scalar */
		s = mnode_build(ST_INTEGER, str);
		if (saml_errno(s) != 0)
			return s;
		mn = create_tensor(0, 1);
	} else {
		/* It has the form literal!index */
		index = atoi(bang+1);
		litlen = bang - str;
		litname = alloca(litlen+1);
		memcpy(litname, str, litlen);
		litname[litlen] = '\0';
		s = mnode_build(ST_LITERAL, litname);
		mn = create_tensor(1, 1);
		range = IRANGE(mn);
		range->lit = s;
		range->start = index;
		range->len = 1;
		s = mnode_build(ST_INTEGER, "1");
	}
	T_DATA(mn)[0] = s;
	return mn;
}

static s_mnode* tensor2tensor (s_mnode *mn1, s_mnode *model)
{
	int nbind = INDICES(mn1), nb_elts;
	s_mnode **data, **data1, *mn;

	if (!model)
		return copy_mnode(mn1);

	nb_elts = count_tensor_elements(mn1);
	mn = create_tensor(nbind, nb_elts);
	copy_index_range(mn, mn1);
	data = T_DATA(mn); data1 = T_DATA(mn1);
	while (nb_elts--)
	    *data++ = mnode_promote(*data1++, *T_DATA(model));
	return mn;
}

static s_mnode* tensor_make (s_mnode* scalar)
{
	s_mnode *mn = create_tensor(0, 1);
	SCALAR(mn) = copy_mnode(scalar);
	return mn;
}

static s_mnode* tensor_zero (s_mnode* model)
{
	int i, nbind = INDICES(model);
	s_mnode *mn;
	index_range *range, *orange;

	/* We create another tensor with the same indices, but with only
	 * one coordinate, equal to zero. */
	mn = create_tensor(nbind, 1);
	orange = IRANGE(model);
	range = IRANGE(mn);
	for (i = 0; i < nbind; i++) {
		range[i].lit = copy_mnode(orange[i].lit);
		range[i].start = orange[i].start;
		range[i].len = 1;	/* smash it! */
	}
	T_DATA(mn)[0] = mnode_zero(T_DATA(model)[0]);
	return mn;
}

static s_mnode* tensor_one (s_mnode* model)
{
	s_mnode *mn = create_tensor(0, 1);

	/* Returns the scalar whose value is '1' */
	SCALAR(mn) = mnode_one(T_DATA(model)[0]);
	return mn;
}

static s_mnode* tensor_negate (s_mnode* mn1)
{
	s_mnode *mn, **data, **data1;
	int nb_elts, nbind = INDICES(mn1);

	nb_elts = count_tensor_elements(mn1);
	mn = create_tensor(nbind, nb_elts);
	copy_index_range(mn, mn1);
	data1 = T_DATA(mn1); data = T_DATA(mn);
	while (nb_elts--)
		*data++ = mnode_negate(*data1++);
	return mn;
}

static gr_string* tensor_stringify (s_mnode* mn)
{
	int i, terms = 0, elt, nb_elts, nbind = INDICES(mn), pos, *position;
	gr_string *grs, *grterm;
	index_range *range = IRANGE(mn);
	s_mnode **data;
	char buffer[24];

	/* Is it a scalar? */
	if (nbind == 0)
		return mnode_stringify(T_DATA(mn)[0]);
	grs = new_gr_string(0);
	nb_elts = count_tensor_elements(mn);
	position = alloca(nbind * sizeof(int));
	data = T_DATA(mn);
	for (elt = 0; elt < nb_elts; elt++, data++) {
		pos = elt;
		for (i = nbind-1; i >= 0; i--) {
			position[i] = pos % range[i].len + range[i].start;
			pos /= range[i].len;
		}
		assert(pos == 0);
		if (!mnode_notzero(*data))
			continue;
		/* So we have a non-zero coordinate */
		grs = grs_append(grs, "+(", 2);
		grterm = mnode_stringify(*data);
		grs = grs_append(grs, grterm->s, grterm->len);
		free(grterm);
		grs = grs_append1(grs, ')');
		/* Now issue the literal!index factors */
		for (i = 0; i < nbind; i++) {
			grs = grs_append1(grs, '.');
			grterm = mnode_stringify(range[i].lit);
			grs = grs_append(grs, grterm->s, grterm->len);
			free(grterm);
			grs = grs_append1(grs, '!');
			sprintf(buffer, "%d", position[i]);
			grs = grs_append(grs, buffer, strlen(buffer));
		}
		terms++;
	}
	if (terms == 0) {
		/* It was zero, but not the scalar zero. */
		grs = grs_append1(grs, '(');
		grterm = mnode_stringify(*T_DATA(mn));
		grs = grs_append(grs, grterm->s, grterm->len);
		free(grterm);
		grs = grs_append1(grs, ')');
		for (i = 0; i < nbind; i++) {
			grs = grs_append1(grs, '.');
			grterm = mnode_stringify(range[i].lit);
			grs = grs_append(grs, grterm->s, grterm->len);
			free(grterm);
			grs = grs_append(grs, "!1", 2);
		}
	}
	return grs;
}

static s_mnode* tensor_add (s_mnode* mn1, s_mnode* mn2)
{
	int i, elt, nb_elts = 1, nbind = INDICES(mn1), pos;
	index_range *range, *range1 = IRANGE(mn1), *range2 = IRANGE(mn2);
	int *idstart, *idlen, *position;
	s_mnode *mn, **data, **entry1, **entry2;

	/* There should be the same number of indices */
	if (nbind != INDICES(mn2))
		return mnode_error(SE_HETERO, "tensor_add");
	idstart = alloca(nbind * sizeof(int));
	idlen = alloca(nbind * sizeof(int));
	for (i = 0; i < nbind; i++) {
		int min1, min2, max1, max2;
		if (range1[i].lit != range2[i].lit)
			return mnode_error(SE_HETERO, "tensor_add");
		min1 = range1[i].start;
		max1 = min1 + range1[i].len;
		min2 = range2[i].start;
		max2 = min2 + range2[i].len;
		if (min2 < min1)
			min1 = min2;
		if (max2 > max1)
			max1 = max2;
		idstart[i] = min1;
		nb_elts *= (idlen[i] = max1 - min1);
	}
	/* All indices are the same. So we can proceed. */
	mn = create_tensor(nbind, nb_elts);
	range = IRANGE(mn);
	for (i = 0; i < nbind; i++) {
		range[i].lit = copy_mnode(range1[i].lit);
		range[i].start = idstart[i];
		range[i].len = idlen[i];
	}
	position = alloca(nbind * sizeof(int));
	data = T_DATA(mn);
	for (elt = 0; elt < nb_elts; elt++, data++) {
		pos = elt;
		for (i = nbind-1; i >= 0; i--) {
			int denom = idlen[i], remainder = pos % denom;
			pos /= denom;
			position[i] = remainder + idstart[i];
			/*
			 * Important optimization here: if (remainder!=0)
			 * then the more significant digits, i.e., the
			 * values of position[i] for lesser values of i,
			 * will be the same for elt and elt-1. So it's not
			 * necessary to update them, so we exit the loop and
			 * avoid some expensive divisions. This optimization
			 * is also used in tensor_mul().
			 */
			if (remainder) break;
		}
		entry1 = get_tensor_entry(mn1, position);
		entry2 = get_tensor_entry(mn2, position);
		if (entry1 != NULL) {
			if (entry2 != NULL)
			    *data = mnode_add(*entry1, *entry2);
			else
			    *data = copy_mnode(*entry1);
		} else {
			/* First operand is zero */
			if (entry2 != NULL)
			    *data = copy_mnode(*entry2);
			else
			    *data = mnode_zero(*T_DATA(mn1));
		}
	}
	return mn;
}

static int intersect_intervals (int start1, int len1, int start2, int len2,
	int *pstart, int *plen)
{
	int tmp, length;

	if (start1 > start2) {
		tmp = start1; start1 = start2; start2 = tmp;
		tmp = len1; len1 = len2; len2 = tmp;
	}
	/* Now we have: start1 <= start2 */
	length = start1 + len1 - start2;
	if (length < 0)
		length = 0;	/* disjoint intervals */
	if (len2 < length)
		length = len2;
	*pstart = start2;
	*plen = length;
	return length;
}

static s_mnode* tensor_mul (s_mnode* mn1, s_mnode* mn2)
{
	int nbind1 = INDICES(mn1), nbind2 = INDICES(mn2), i, elt, nb_elts,
		elt_c, nb_elts_c, pos, pos_c, nb_f_ind = 0, nb_c_ind = 0;
	index_range *range1 = IRANGE(mn1), *range2 = IRANGE(mn2), *p1, *p2,
		*range_c, *range, *range1_end, *range2_end, *pc, *pf;
	s_mnode *prod, *zero, **data, *entry1, *entry2, *pt1, *pt2;
	s_mnode **base_entry1, **base_entry2;
	int *position1, *position2;
	int **origpos, **origpos1, **origpos2;
	int *offset1, *offset2;

#if 0
	/* Test if one of the factors is a scalar */
	if (INDICES(mn1) == 0)
		return tensor_mul_scalar(mn2, SCALAR(mn1), 0);
	else if (INDICES(mn2) == 0)
		return tensor_mul_scalar(mn1, SCALAR(mn2), 0);
#endif
	/* First step: find the common indices between mn1 and mn2 */
	pf = range = alloca((nbind1+nbind2) * sizeof(index_range));
	origpos = alloca((nbind1+nbind2) * sizeof(int*));
	i = (nbind1 < nbind2) ? nbind1 : nbind2;
	pc = range_c = alloca(i * sizeof(index_range));
	origpos1 = alloca(i * sizeof(int*));
	origpos2 = alloca(i * sizeof(int*));
	position1 = alloca(nbind1 * sizeof(int));
	position2 = alloca(nbind2 * sizeof(int));
	range1_end = range1 + nbind1;
	range2_end = range2 + nbind2;
	p1 = range1; p2 = range2; nb_elts = nb_elts_c = 1;
	/* We assume the literals are in ascending order */
	while (p1 < range1_end && p2 < range2_end) {
		s_mnode *lit1 = p1->lit, *lit2 = p2->lit;
		if (lit1 == lit2) {
			/* Common literal found */
			pc->lit = lit1;
			if (!intersect_intervals(p1->start, p1->len,
			  p2->start, p2->len, &(pc->start), &(pc->len)))
				nb_elts = 0;	/* No intersection */
			nb_elts_c *= pc->len;
			origpos1[pc-range_c] = position1 + (p1-range1);
			origpos2[pc-range_c] = position2 + (p2-range2);
			++pc, ++p1, ++p2;
		} else if (lit1 < lit2) {
			/* The literal in mn1 doesn't appear in mn2 */
			pf->lit = copy_mnode(lit1);
			pf->start = p1->start;
			nb_elts *= (pf->len = p1->len);
			origpos[pf-range] = position1 + (p1-range1);
			++pf, ++p1;
		} else {
			/* The literal in mn2 doesn't appear in mn1 */
			pf->lit = copy_mnode(lit2);
			pf->start = p2->start;
			nb_elts *= (pf->len = p2->len);
			origpos[pf-range] = position2 + (p2-range2);
			++pf, ++p2;
		}
	}
	while (p1 < range1_end) {
		pf->lit = copy_mnode(p1->lit);
		pf->start = p1->start;
		nb_elts *= (pf->len = p1->len);
		origpos[pf-range] = position1 + (p1-range1);
		++pf, ++p1;
	}
	while (p2 < range2_end) {
		pf->lit = copy_mnode(p2->lit);
		pf->start = p2->start;
		nb_elts *= (pf->len = p2->len);
		origpos[pf-range] = position2 + (p2-range2);
		++pf, ++p2;
	}
	nb_c_ind = pc - range_c;	/* Indices to contract */
	nb_f_ind = pf - range;		/* Other indices */
#if 0
	fprintf(stderr, "tensor_mul: %d free indices, %d contracted\n",
		nb_f_ind, nb_c_ind);
#endif
	assert(nb_c_ind*2 + nb_f_ind == nbind1 + nbind2);
	prod = create_tensor(nb_f_ind, (nb_elts? nb_elts : 1));
	memcpy(IRANGE(prod), range, nb_f_ind * sizeof(index_range));
	zero = mnode_zero(T_DATA(mn1)[0]);
	if (nb_elts == 0) {
		/* For some contracted index, the ranges didn't overlap,
		   thus the result is zero. */
		range = IRANGE(prod);
		for (i = 0; i < nb_f_ind; i++)
			range[i].len = 1;
		T_DATA(prod)[0] = zero;
		return prod;
	}
	/* Compute the offsets when the contracted indices move */
	offset1 = alloca(nb_elts_c * sizeof(int));
	offset2 = alloca(nb_elts_c * sizeof(int));
	memset(position1, 0, nbind1 * sizeof(int));
	memset(position2, 0, nbind2 * sizeof(int));
	for (elt_c = 0; elt_c < nb_elts_c; elt_c++) {
	  pos_c = elt_c;
	  for (i = nb_c_ind-1; i >= 0; i--) {
	    *(origpos1[i]) = *(origpos2[i]) = pos_c % range_c[i].len;
	    pos_c /= range_c[i].len;
	  }
	  pos = 0;
	  for (i = 0; i < nbind1; i++)
	    pos = pos * range1[i].len + position1[i];
	  offset1[elt_c] = pos;
	  pos = 0;
	  for (i = 0; i < nbind2; i++)
	    pos = pos * range2[i].len + position2[i];
	  offset2[elt_c] = pos;
	}
	/* Now compute all individual elements */
	zero->refs += (nb_elts - 1);
	data = T_DATA(prod);
	for (elt = 0; elt < nb_elts; elt++,data++) {
	  pos = elt;
	  for (i = nb_f_ind-1; i >= 0; i--) {
	    int denom = range[i].len, remainder = pos % denom;
	    pos /= denom;
	    *(origpos[i]) = remainder + range[i].start;
	    if (remainder)
	    	break;	   /* same optimization as in tensor_add() */
	  }
	  for (i = 0; i < nb_c_ind; i++)
	    *(origpos1[i]) = *(origpos2[i]) = range_c[i].start;
	  base_entry1 = get_tensor_entry(mn1, position1);
	  base_entry2 = get_tensor_entry(mn2, position2);
	  *data = zero;
	  for (elt_c = 0; elt_c < nb_elts_c; elt_c++) {
	    entry1 = base_entry1[offset1[elt_c]];
	    if (!mnode_notzero(entry1))
	    	continue;
	    entry2 = base_entry2[offset2[elt_c]];
	    if (!mnode_notzero(entry2))
	    	continue;
	    pt1 = mnode_mul(entry1, entry2);
	    pt2 = mnode_add(*data, pt1);
	    unlink_mnode(*data); unlink_mnode(pt1);
	    *data = pt2;
	  }
	}
	return prod;
}

static s_mnode* tensor_div (s_mnode* mn1, s_mnode* mn2)
{
	/* The denominator must be a scalar */
	if (INDICES(mn2))
		return mnode_error(SE_NSCALAR, "tensor_div");
	return tensor_mul_scalar(mn1, SCALAR(mn2), 1);
}

static s_mnode* tensor_gcd (s_mnode* mn1, s_mnode *mn2)
{
	s_mnode *s, *t;
	
	/* Both operands must be scalars */
	if (INDICES(mn1) || INDICES(mn2))
		return mnode_error(SE_NSCALAR, "tensor_gcd");
	s = mnode_gcd(SCALAR(mn1), SCALAR(mn2));
	t = tensor_make(s);
	unlink_mnode(s); return t;
}

static int tensor_notzero (s_mnode* mn)
{
	int nb_elt = count_tensor_elements(mn);
	s_mnode **data = T_DATA(mn);

	while (nb_elt--)
	    if (mnode_notzero(*data++))
	    	return 1;
	/* All coordinates are equal to zero */
	return 0;
}

s_mnode* tensor_move_literal (s_mnode* mn1, s_mnode* l1, s_mnode* l2)
{
	int i, elt, nb_elts, nbind, nbind1 = INDICES(mn1), pos, pos1, pos2,
		c_start, c_len;
	index_range *range1 = IRANGE(mn1), *p, *p1, *range;
	s_mnode *mn, **data, **entry, *t1, *zero;
	int *origpos, *position;

	pos1 = pos2 = -1;
	for (i = 0; i < nbind1; i++) {
		if (range1[i].lit == l1)
			pos1 = i;
		else if (range1[i].lit == l2)
			pos2 = i;
	}
	/* The first literal MUST appear in the tensor */
	if (pos1 < 0)
		return mnode_error(SE_INDEX, "tensor_move_literal");
	if (pos2 < 0) {
		/*
		 * The second literal doesn't appear in mn1, so it's really
		 * only a renaming. The only problem is that l2 may not have
		 * the same position as l1 WRT the other literals, which will
		 * force us to some data shuffling.
		 */
		nbind = nbind1;
		nb_elts = count_tensor_elements(mn1);
		mn = create_tensor(nbind, nb_elts);
		p = range = IRANGE(mn);
		origpos = alloca(nbind * sizeof(int));
		position = alloca(nbind1 * sizeof(int));
		for (p1 = range1; p1 < range1+nbind1; p1++) {
			if (p1 == range1+pos1)
				continue;
			if (l2 && p1->lit > l2) {
				/* Insert the second literal here */
				p->lit = copy_mnode(l2);
				p->start = range1[pos1].start;
				p->len = range1[pos1].len;
				origpos[p-range] = pos1;
				p++;
				/* Nuke l2, so it won't be inserted again */
				l2 = NULL;
			}
			p->lit = copy_mnode(p1->lit);
			p->start = p1->start;
			p->len = p1->len;
			origpos[p-range] = p1-range1;
			p++;
		}
		if (p-range != nbind) {
			assert(p-range == nbind-1 && l2 != NULL);
			p->lit = copy_mnode(l2);
			p->start = range1[pos1].start;
			p->len = range1[pos1].len;
			origpos[p-range] = pos1;
		}
		/* Good. Now copy the data */
		data = T_DATA(mn);
		for (elt = 0; elt < nb_elts; elt++,data++) {
			pos = elt;
			for (i = nbind-1; i >= 0; i--) {
				position[origpos[i]] = pos % range[i].len
					+ range[i].start;
				pos /= range[i].len;
			}
			assert(pos == 0);
			entry = get_tensor_entry(mn1, position);
			assert(entry != NULL);
			*data = copy_mnode(*entry);
		}
		return mn;
	}
	/* Second case: both indices appear, so we must contract them. */
	nbind = nbind1 - 2;
	origpos = alloca(nbind * sizeof(int));
	position = alloca(nbind1 * sizeof(int));
	nb_elts = 1;
	for (i = pos = 0; i < nbind1; i++)
		if (i != pos1 && i != pos2) {
			origpos[pos++] = i;
			nb_elts *= range1[i].len;
		}
	assert(pos == nbind);
	if (!intersect_intervals(range1[pos1].start, range1[pos1].len,
	    range1[pos2].start, range1[pos2].len, &c_start, &c_len))
	    	nb_elts = 0;
	mn = create_tensor(nbind, (nb_elts ? nb_elts : 1));
	range = IRANGE(mn);
	for (i = 0; i < nbind; i++) {
		int op = origpos[i];
		range[i].lit = copy_mnode(range1[op].lit);
		range[i].start = range1[op].start;
		range[i].len = (nb_elts ? range1[op].len : 1);
	}
	zero = mnode_zero(*T_DATA(mn1));
	if (nb_elts == 0) {
		*T_DATA(mn) = zero;
		return mn;
	}
	zero->refs += nb_elts - 1;
	data = T_DATA(mn);
	for (elt = 0; elt < nb_elts; elt++,data++) {
		pos = elt;
		for (i = nbind-1; i >= 0; i--) {
			position[origpos[i]] = pos % range[i].len
				+ range[i].start;
			pos /= range[i].len;
		}
		assert(pos == 0);
		*data = zero;
		/* Now iterate over the contracted indices */
		for (i = c_start; i < c_start+c_len; i++) {
			position[pos1] = position[pos2] = i;
			entry = get_tensor_entry(mn1, position);
			if (entry == NULL)
				continue;
			t1 = mnode_add(*data, *entry);
			unlink_mnode(*data);
			*data = t1;
		}
	}
	return mn;
}

s_mnode* tensor_subs (s_mnode* mn1, s_mnode* e1, s_mnode* e2)
{
	s_mnode *mn, **data, **data1;
	int nbind = INDICES(mn1), nb_elts;

	if (e1->type == ST_TENSOR) {
		if (INDICES(e1))
			return mnode_error(SE_NSCALAR, "tensor_subs");
		e1 = SCALAR(e1);
	}
	if (e2->type == ST_TENSOR) {
		if (INDICES(e2))
			return mnode_error(SE_NSCALAR, "tensor_subs");
		e2 = SCALAR(e2);
	}
	nb_elts = count_tensor_elements(mn1);
	mn = create_tensor(nbind, nb_elts);
	copy_index_range(mn, mn1);
	data = T_DATA(mn); data1 = T_DATA(mn1);
	while (nb_elts--)
		*data++ = mnode_subs(*data1++, e1, e2);
	return mn;
}

s_mnode* tensor_diff (s_mnode* mn1, s_mnode* lit)
{
	s_mnode *mn, **data1, **data;
	int nbind = INDICES(mn1), nb_elts;
	
	if (lit->type == ST_TENSOR) {
		if (INDICES(lit))
			return mnode_error(SE_NSCALAR, "tensor_diff");
		lit = SCALAR(lit);
	}
	nb_elts = count_tensor_elements(mn1);
	mn = create_tensor(nbind, nb_elts);
	copy_index_range(mn, mn1);
	data = T_DATA(mn); data1 = T_DATA(mn1);
	while (nb_elts--)
		*data++ = mnode_diff(*data1++, lit);
	return mn;
}

s_mnode* tensor_sylvester (s_mnode* t1, s_mnode* t2, s_mnode *lit)
{
	s_mnode *s, *t;

	if (INDICES(t1) || INDICES(t2))
		return mnode_error(SE_NSCALAR, "tensor_sylvester");
	t1 = SCALAR(t1);
	t2 = SCALAR(t2);
	if (lit->type == ST_TENSOR && !INDICES(lit))
		lit = SCALAR(lit);
	s = mnode_elim(lit, t1, t2);
	t = tensor_make(s);
	unlink_mnode(s);
	return t;
}