/*
 * Copyright 1995,96 Thierry Bousch
 * Licensed under the Gnu Public License, Version 2
 *
 * $Id: mnode.c,v 3.6 1996/08/18 09:34:08 bousch Exp $
 *
 * Generic operations on mnodes; they are dispatched to the real operations,
 * according to the dynamic types of the mnodes.
 */

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

long nb_mnodes_allocated = 0;
long nb_mnodes_freed = 0;
long nb_mnodes_reserved = 0;

s_mtype* mtype_table[MAX_MTYPES];

void register_mtype (int t, void* desc)
{
	char msg[100];

	if (!MTYPE_OK(t)) {
		sprintf(msg, "math type %d out of bounds", t);
		saml_panic(msg);
	}
	if (MTYPE_INUSE(t)) {
		sprintf(msg, "math type %d already registered", t);
		saml_panic(msg);
	}
#ifdef DEBUG_MTYPES
	if (getenv("DEBUG_MTYPES"))
		fprintf(stderr, "Adding type %3d: \"%s\"\n", t,
		((s_mtype*)desc)->name);
#endif
	mtype_table[t] = desc;
}

/*
 * Unary operations
 */

#define RETURN_OP(where,op,args)			\
	return (op)? (*op)args : mnode_error(SE_ONSUPP,where)

#define UNARY_OPERATION(opname)				\
  s_mnode* mnode_##opname (s_mnode *n)			\
  {							\
	s_mnode* (*op)(s_mnode*);			\
	op = mtype_table[n->type]->opname;		\
	RETURN_OP(#opname, op, (n));			\
  }

UNARY_OPERATION(zero)
UNARY_OPERATION(negate)
UNARY_OPERATION(one)
UNARY_OPERATION(invert)
UNARY_OPERATION(sqrt)

/*
 * Binary operations
 */

#define BINARY_OPERATION(opname)				\
  s_mnode* mnode_##opname (s_mnode *n1, s_mnode *n2)		\
  {								\
  	s_mnode* (*op)(s_mnode*, s_mnode*);			\
  	int t1 = n1->type;					\
	int t2 = n2->type;					\
								\
	if (t1 == t2) {						\
		op = mtype_table[t1]->opname;			\
		RETURN_OP(#opname, op, (n1,n2));		\
	}							\
	else if (t1 == ST_VOID)					\
		return copy_mnode(n1);				\
	else if (t2 == ST_VOID)					\
		return copy_mnode(n2);				\
	/* Type conflict */					\
	return mnode_error(SE_TCONFL, #opname);			\
  }

BINARY_OPERATION(add)
BINARY_OPERATION(sub)
BINARY_OPERATION(mul)
BINARY_OPERATION(div)
BINARY_OPERATION(gcd)

/*
 * Boolean functions; return -1 if the operation is undefined
 */

int mnode_notzero (s_mnode *n)
{
	int (*op)(s_mnode*);
	op = mtype_table[n->type]->notzero;
	return op? (*op)(n) : -1;
}

int mnode_isneg (s_mnode *n)
{
	int (*op)(s_mnode*);
	op = mtype_table[n->type]->isneg;
	return op? (*op)(n) : -1;
}

int mnode_info (s_mnode *n, int param)
{
	int (*op)(s_mnode*, int);
	op = mtype_table[n->type]->info;
	return op? (*op)(n,param) : -1;
}

int mnode_differ (s_mnode *n1, s_mnode *n2)
{
	int (*op)(s_mnode*, s_mnode*);
	int t1 = n1->type, t2 = n2->type;

	if (t1 == t2) {
		op = mtype_table[t1]->differ;
		if (op)
		    return (*op)(n1,n2);
	}
	return -1;
}

int mnode_lessthan (s_mnode *n1, s_mnode *n2)
{
	int (*op)(s_mnode*, s_mnode*);
	int t1 = n1->type, t2 = n2->type;

	if (t1 == t2) {
		op = mtype_table[t1]->lessthan;
		if (op)
		    return (*op)(n1,n2);
	}
	return -1;
}

/*
 * Conversion from/to a string
 */

s_mnode* mnode_build (int typeid, const char* string)
{
	s_mnode* (*op)(const char*);

	if (!string)
		return mnode_error(SE_STRING, "build");
	if (!MTYPE_OK(typeid) || !MTYPE_INUSE(typeid))
		return mnode_error(SE_NSTYPE, "build");
	op = mtype_table[typeid]->build;
	RETURN_OP("build", op, (string));
}

s_mnode* mnode_make (int typeid, s_mnode* n)
{
	s_mnode* (*op)(s_mnode*);

	if (!MTYPE_OK(typeid) || !MTYPE_INUSE(typeid))
		return mnode_error(SE_NSTYPE, "make");
	op = mtype_table[typeid]->make;
	RETURN_OP("make", op, (n));
}

s_mnode* mnode_etc (s_mnode *n, int param, void *data)
{
	s_mnode* (*op)(s_mnode*, int, void*);
	op = mtype_table[n->type]->etc;
	RETURN_OP("etc", op, (n,param,data));
}

gr_string* mnode_stringify (s_mnode *n)
{
	gr_string* (*op)(s_mnode*);
	op = mtype_table[n->type]->stringify;
	return op? (*op)(n) : NULL;
}

/*
 * Some very standard operations
 */

s_mnode* mnode_mod (s_mnode *n1, s_mnode *n2)
{
	s_mnode *quot, *tmp, *rem;

	quot = mnode_div(n1, n2);
	if (quot->type == ST_VOID) {
		/* If we can't divide, just return n1. */
		unlink_mnode(quot);
		return copy_mnode(n1);
	}
	tmp = mnode_mul(quot, n2);
	unlink_mnode(quot);
	rem = mnode_sub(n1, tmp);
	unlink_mnode(tmp);
	return rem;
}

s_mnode* mnode_power (s_mnode *n, int expo)
{
	s_mnode *n1, *n2, *n3;
	extern int apoly_length(s_mnode*);

	if (!expo)
		return mnode_one(n);
	if (expo < 0) {
		n1 = mnode_invert(n);
		if (n1->type == ST_VOID)
			return n1;
		n2 = mnode_power(n1, -expo);
		unlink_mnode(n1);
		return n2;
	}
	/*
	 * Here we know that (expo > 0). If the operand is a polynomial,
	 * it's quicker to do multiplications one-by-one, because
	 * multiplication of long polynomials is very time-consuming.
	 * In all other cases, we use the binary method.
	 */
	if (expo < 4 ||
	    (n->type == ST_POLY  && ((smn_ptr)n)->length > 2) ||
	    (n->type == ST_APOLY &&    apoly_length(n)   > 1) ||
	    (n->type == ST_UPOLY && ((smn_ptr)n)->length > 1)) {
		n1 = copy_mnode(n);
		while (--expo) {
			n2 = mnode_mul(n1, n);
			unlink_mnode(n1);
			n1 = n2;
		}
		return n1;
	}
	/* Binary method */
	n1 = mnode_one(n);
	n2 = copy_mnode(n);
	while (1) {
		if (expo & 1) {
			n3 = mnode_mul(n1, n2);
			unlink_mnode(n1);
			n1 = n3;
		}
		if ((expo >>= 1) == 0)
			break;
		n3 = mnode_mul(n2, n2);
		unlink_mnode(n2);
		n2 = n3;
	}
	unlink_mnode(n2);
	return n1;
}

s_mnode* mn_std_sub (s_mnode *n1, s_mnode *n2)
{
	/* Standard substraction in a commutative group */
	s_mnode* n3 = mnode_negate(n2);
	s_mnode* n4 = mnode_add(n1,n3);
	unlink_mnode(n3);
	return n4;
}

s_mnode* mn_std_div (s_mnode *n1, s_mnode *n2)
{
	/* Standard division in a commutative field */
	s_mnode* n3 = mnode_invert(n2);
	s_mnode* n4 = mnode_mul(n1,n3);
	unlink_mnode(n3);
	return n4;
}

int mn_std_differ (s_mnode *n1, s_mnode *n2)
{
	/* Standard comparison in a group */
	s_mnode* n3 = mnode_sub(n1,n2);
	int s = mnode_notzero(n3);
	unlink_mnode(n3); return s;
}

int mn_std_lessthan (s_mnode *n1, s_mnode *n2)
{
	s_mnode* n3 = mnode_sub(n1,n2);
	int s = mnode_isneg(n3);
	unlink_mnode(n3); return s;
}

s_mnode* mn_field_gcd (s_mnode *n1, s_mnode *n2)
{
	if (mnode_notzero(n1) || mnode_notzero(n2))
		return mnode_one(n1);
	else return mnode_zero(n1);
}

s_mnode* mn_euclidean_gcd (s_mnode *n1, s_mnode *n2)
{
	s_mnode *arg1, *arg2, *quot, *tmp, *rem;

	/* An important optimization */
	if (n1 == n2 || !mnode_differ(n1,n2))
		return copy_mnode(n1);

	arg1 = copy_mnode(n1);
	arg2 = copy_mnode(n2);
	while (mnode_notzero(arg2)) {
		quot = mnode_div(arg1, arg2);
		if (mnode_notzero(quot)) {
			tmp = mnode_mul(arg2, quot);
			rem = mnode_sub(arg1, tmp);
			unlink_mnode(tmp);
			unlink_mnode(arg1);
		} else {
			rem = arg1;
		}
		unlink_mnode(quot);
		arg1 = arg2;
		arg2 = rem;
	}
	unlink_mnode(arg2);
	return arg1;
}

/*
 * Miscillaneous functions
 */

void destroy_mnode (s_mnode *n)
{
	(mtype_table[n->type]->free)(n);
	++nb_mnodes_freed;
}

void mstd_free (std_mnode *n)
{
	int len;
	s_mnode **p;

	for (len = n->length, p = n->x; len; ++p, --len)
		unlink_mnode(*p);
	free(n);
}

s_mnode* mnode_det (s_mnode* matrix)
{
	extern s_mnode* matrix_determinant(s_mnode*);

	if (matrix->type == ST_MATRIX)
		return matrix_determinant(matrix);
	return mnode_error(SE_ONSUPP, "determinant");
}

s_mnode* mnode_diff (s_mnode* expr, s_mnode* lit)
{
	extern s_mnode* poly_diff(s_mnode*, s_mnode*);
	extern s_mnode* apoly_diff(s_mnode*, s_mnode*);
	extern s_mnode* tensor_diff(s_mnode*, s_mnode*);

	switch (expr->type) {
	  case ST_POLY:
		return poly_diff(expr, lit);
	  case ST_APOLY:
	  	return apoly_diff(expr, lit);
	  case ST_TENSOR:
		return tensor_diff(expr, lit);
	}
	return mnode_error(SE_ONSUPP, "diff");
}

s_mnode* mnode_move_lit (s_mnode* tensor, s_mnode* lit1, s_mnode* lit2)
{
	extern s_mnode* tensor_move_literal(s_mnode*, s_mnode*, s_mnode*);

	if (tensor->type == ST_TENSOR) {
		if (lit1->type == ST_LITERAL && lit2->type == ST_LITERAL)
			return tensor_move_literal(tensor, lit1, lit2);
		return mnode_error(SE_TCONFL, "move_lit");
	}
	return mnode_error(SE_ONSUPP, "move_lit");
}

s_mnode* mnode_subs (s_mnode* expr, s_mnode* e1, s_mnode* e2)
{
	extern s_mnode* poly_subs(s_mnode*, s_mnode*, s_mnode*);
	extern s_mnode* apoly_subs(s_mnode*, s_mnode*, s_mnode*);
	extern s_mnode* tensor_subs(s_mnode*, s_mnode*, s_mnode*);

	switch(expr->type) {
	  case ST_POLY:
		return poly_subs(expr, e1, e2);
	  case ST_APOLY:
	  	return apoly_subs(expr, e1, e2);
	  case ST_TENSOR:
		return tensor_subs(expr, e1, e2);
	}
	return mnode_error(SE_ONSUPP, "subs");
}

s_mnode* mnode_elim (s_mnode* lit, s_mnode* e1, s_mnode* e2)
{
	extern s_mnode* poly_sylvester(s_mnode*, s_mnode*, s_mnode*);
	extern s_mnode* apoly_sylvester(s_mnode*, s_mnode*, s_mnode*);
	extern s_mnode* tensor_sylvester(s_mnode*, s_mnode*, s_mnode*);

	if (e1->type != e2->type)
		return mnode_error(SE_TCONFL, "elim");
	switch(e1->type) {
	    case ST_POLY:
		return poly_sylvester(e1, e2, lit);
	    case ST_APOLY:
		return apoly_sylvester(e1, e2, lit);
	    case ST_TENSOR:
		return tensor_sylvester(e1, e2, lit);
	}
	return mnode_error(SE_ONSUPP, "elim");
}