/*  This file handles buit in operators & functions.
    Copyright (C) 1990  Marty White

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License
    as published by the Free Software Foundation; either version 2
    of the License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

#include <stdio.h>
#include <math.h>
#include "compiler.h"

#include "physcalc.h"
#include "physdecl.h"

#ifdef __PROTOTYPES__
LOCAL void change_type(NODEP n, int type);
LOCAL NODEP derivative(NODEP n);
#endif

/*----- Miscelaneous Math related routines -----*/

EXPORT int number(n)	/* is n a number? */
NODEP n;
{
	return n->type!=SNODE && n->type<LNODE;
}

LOCAL void change_type(n,type)	/* Convert node N to type TYPE */
NODEP n;		/* Assumes conversion is 'upward', but will also */
int type;		/* force conversion from any type to a double */
{
	int i;
	double x;

	switch (n->type) {
	case INODE:
		switch (type) {
		case FNODE:
			i = (int) n->data->lng;
			reallocnode(n,FNODE,sizeof(struct fractstruct));
			n->data->frt.numer = i;
			n->data->frt.denom = 1;
			break;
		case DNODE:
			x = (double) n->data->lng;
			reallocnode(n,DNODE,sizeof(double));
			n->data->dbl = x;
			break;
		case NNODE:
			x = (double) n->data->lng;
			reallocnode(n,NNODE,sizeof(DNUM));
			erasednum(&n->data->dmn);
			n->data->dmn.num = x;
			break;
		}
		break;
	case FNODE:
		x = (double) n->data->frt.numer / (double) n->data->frt.denom;
		if (type==DNODE) {
			reallocnode(n,DNODE,sizeof(double));
			n->data->dbl = x;
		} else {
			reallocnode(n,DNODE,sizeof(DNUM));
			erasednum(&n->data->dmn);
			n->data->dmn.num = x;
		}
		break;
	case DNODE:
		if (type==NNODE) {
			x = n->data->dbl;
			reallocnode(n,NNODE,sizeof(DNUM));
			erasednum(&n->data->dmn);
			n->data->dmn.num = x;
		}
		break;
	case NNODE:
		if (type==DNODE) {
			x = n->data->dmn.num;
			reallocnode(n,DNODE,sizeof(double));
			n->data->dbl = x;
		}
	}
}

EXPORT void force_same_type(x,y)	/* Change x & y to be the same type */
NODEP x,y;
{
	int rtype=INODE;

	/* Determine common type (assume integer) */

	if (x->type==FNODE || y->type==FNODE)
		rtype=FNODE;
	if (x->type==DNODE || y->type==DNODE)
		rtype=DNODE;
	if (x->type==NNODE || y->type==NNODE)
		rtype=NNODE;

	/* Now convert each to rtype */
	if (x->type != rtype)
		change_type(x,rtype);
	if (y->type != rtype)
		change_type(y,rtype);
}

EXPORT void simplify_frt(n)	/* Simplify a fraction */
NODEP n;
{
	int a,b,r;

	if (a = n->data->frt.numer) {
		b = n->data->frt.denom;
		if (a<0)
			a *= -1;
		if (a<b) {					/* Find g.c.d. using Euclid's algorithm */
			r = a;		/* make (a) the larger */
			a = b;
			b = r;
		}
		do {
			r = a % b;	/* remainder after a/b */
			a = b;
			b = r;
		} while (r);	/* until 0 remainder */
		if (a>1) {
			n->data->frt.numer /= a;
			n->data->frt.denom /= a;
		}
	} else
		n->data->frt.denom = 1;
	if (n->data->frt.denom==1) {	/* Change to int if denominator is 1 */
		a = n->data->frt.numer;
		reallocnode(n,INODE,sizeof(long));
		n->data->lng = (long) a;
	}
}

EXPORT void factor(s)	/* Output all integer factors of an expression */
char const *s;
{
	NODEP n;
	int i,root,prime;
	long a;

	n = solve(&s);
	if (error) {
		printf("To use FACTOR type:\n  FACTOR <integer value>\n");
		return;
	}
	if (n->type != INODE) {
		printf("Non-integers are un-factorable.\n");
		deallocnode(n);
		return;
	}
	a = n->data->lng;
	deallocnode(n);
	if (!a) {
		printf("The factors of zero are zero and anything!\n");
		return;
	}
	printf("Factors are: 1,%d,",a);
	prime = TRUE;
	root = (int) sqrt((double) a);
	for (i=2; i<=root; i++)
		if ( (a/i)*i == a) {
			printf("%d,%d,",i,a/i);
			prime = FALSE;
		}
	printf("\b \n");
	if (prime)
		printf("This is a prime number.\n");
	else {
		printf("The prime factorization is: ");
		for (i=2; i<=root; i++)
			if ( (a/i)*i == a) {
				printf("%d,",i);
				a = a/i;
				i=1;
				root = sqrt((double) a);
			}
		printf("%d\n",a);
	}
}

/*----- Built in functions -----*/
/* Caller is responsible for de-allocation of parameters */

LOCAL NODEP factorial(x)		/* factorial operator/function */
NODEP x;
{
	int i;
	long n, old_n;

	if (x->type != INODE) {
		error = TYPERR;
		return NULL;
	}
	n = 1;
	old_n = 1;
	for (i=2; i <= x->data->lng; i++) {
		n *= i;
		if (n <= old_n) {
			error = MTHERR;
			return NULL;
		}
		old_n = n;
	}
	return allocinode(n);
}

LOCAL NODEP integrate(f,a,b,p)
	/* integrate f from a to b using trapezoidal rule */
NODEP f,a,b,p;	/* f should be a string; p should be a positive int or NULL */
{
	double deltax,a1,b1,x,y;
	NODEP n,n2,n3,na,nb,np;
	static double parts = 20.0;

	if (!f || !a || !b)
		return NULL;
	if (p) {
		np = evaluate(p);
		if (np->type != INODE || np->data->lng < 1L || np->data->lng>32767L) {
			printf("Partitions must be an integer from 1 to 32768.\n");
			deallocnode(np);
			return NULL;
		}
		parts = (double) np->data->lng;
		deallocnode(np);
	}
	na = evaluate(a);
	if (error) {
		deallocnode(na);
		return NULL;
	}
	nb = evaluate(b);
	if (error) {
		deallocnode(na);
		deallocnode(nb);
		return NULL;
	}
	change_type(na,DNODE);
	change_type(nb,DNODE);
	a1 = na->data->dbl;
	b1 = nb->data->dbl;
	deallocnode(na);
	deallocnode(nb);

	deltax = (b1-a1) / parts;
	printf("Integrating...\n Delta = %lg\n",deltax);
	if (fabs(deltax) == 0.0) {
		printf("Cannot integrate on a non-existent interval!\n");
		return NULL;
	}
	if (parts>30)
		printf(" (Please be patient while I think...)\n");

	n = alloclnode(LNODE);	/* Create a node calling f with a */
	linknode(n, copynode(f) );

	n2 = allocdnode(a1);
	linknode(n,n2);

	n3 = evaluate(n);
	change_type(n3,DNODE);
	y = n3->data->dbl;

	for (x = a1 + deltax; x <= b1 - deltax; x += deltax) {
		n2->data->dbl = x;
		deallocnode(n3);
		n3 = evaluate(n);
		change_type(n3,DNODE);
		y += 2 * n3->data->dbl;
	}

	n2->data->dbl = b1;
	deallocnode(n3);
	n3 = evaluate(n);
	change_type(n3,DNODE);
	y += n3->data->dbl;

	y *= deltax * 0.5;

	deallocnode(n);
	deallocnode(n3);

	return allocdnode(y);
}

LOCAL NODEP derivative(n)	/* derivative function (not too useful) */
NODEP n;
{
	NODEP n2,n3,n4,n5;

	switch (n->type) {
	case SNODE:						/* Derivative of a variable is 1 */
		n2 = allocinode(1L);
		break;
    case 0 /*NULL*/:
		n2 = NULL;
		break;
	case INODE:
	case FNODE:
	case DNODE:
	case NNODE:              		/* derivative of a constant is 0 */
		n2 = allocinode(0L);
		break;
	case LNODE:
		n2 = derivative(n->data->list[1]);
		break;
	default:								/* Implement other rules */
		switch (n->type - LNODE) {
		case 1:							/* Addition */
		case 2:							/* Subtracton */
			n2 = alloclnode2(n->type, derivative(n->data->list[0]), derivative(n->data->list[1]));
			break;
		case 3:							/* multiplication */
			n2 = alloclnode(LNODE + 1);	/* create addition node */
			n3 = alloclnode(LNODE + 3);	/* mult node */
			linknode(n2,n3);
			linknode(n3, derivative(n->data->list[0]) );
			linknode(n3, copynode(n->data->list[1]) );
			n3 = alloclnode(LNODE + 3);
			linknode(n2,n3);
			linknode(n3, copynode(n->data->list[0]) );
			linknode(n3, derivative(n->data->list[1]) );
			break;
		case 4:							/* division */
			n2 = alloclnode(LNODE + 4);
			n3 = alloclnode(LNODE + 2);
			linknode(n2,n3);
			n4 = alloclnode(LNODE + 3);
			linknode(n3,n4);
			linknode(n4, derivative(n->data->list[0]) );
			linknode(n4, copynode(n->data->list[1]) );
			n4 = alloclnode(LNODE + 3);
			linknode(n3,n4);
			linknode(n4, copynode(n->data->list[0]) );
			linknode(n4, derivative(n->data->list[1]) );
			n3 = alloclnode(LNODE + 5);
			linknode(n2,n3);
			linknode(n3, copynode(n->data->list[1]) );
			linknode(n3, allocinode(2L) );
			break;
		case 5:							/* power rule */
			n2 = alloclnode(LNODE + 3);
			n3 = alloclnode(LNODE + 3);	/* deriv. of the outside */
			linknode(n2,n3);
			linknode(n3, copynode(n->data->list[1]) );
			n4 = alloclnode(LNODE + 5);
			linknode(n3,n4);
			linknode(n4, copynode(n->data->list[0]) );
			n5 = alloclnode(LNODE + 2);
			linknode(n4,n5);
			linknode(n5, copynode(n->data->list[1]) );
			linknode(n5, allocinode(1L) );
				/* times the deriv. of the inside */
			linknode(n2, derivative(n->data->list[0]) );
			break;
		default:
			error = TYPERR;
			n2 = NULL;
			break;
		}
		break;
	}
	n3 = evaluate(n2);
	deallocnode(n2);
	return n3;
}


#ifndef __PROTOTYPES__
EXPORT double (*do_funct[FUNCS])()
#else
EXPORT double (*do_funct[FUNCS])(double parm)
#endif
  = {
		exp,log,sqrt,sin,cos,tan,asin,acos,atan,log10
	};


EXPORT char const functions[FUNCS][FNLEN] = {
		"EXP","LN","SQRT","SIN","COS","TAN","ARCSIN","ARCCOS","ARCTAN","LOG"
	};

EXPORT NODEP (*do_funct2[FUNCS2])() = {
		factorial,integrate,derivative
	};

EXPORT char const functions2[FUNCS2][FNLEN] = {
		"FACT","INTEGRATE","DERIVATIVE"
	};

/*---- Arithmetic operator functions ----*/

LOCAL NODEP add(x,y)
NODEP x,y;
{
	NODEP z;

	switch (x->type) {
	case INODE:
		z = allocinode(x->data->lng + y->data->lng);
		break;
	case FNODE:
		z = allocfnode(x->data->frt.numer * y->data->frt.denom +
						y->data->frt.numer * x->data->frt.denom,
						x->data->frt.denom * y->data->frt.denom);
		break;
	case DNODE:
		z = allocdnode(x->data->dbl + y->data->dbl);
		break;
	case NNODE:
		z = allocnnode();
		if (!check_equ_dim(&x->data->dmn,&y->data->dmn))
			printf("Adding unlike dimensions!\n");
		else
			copydnum(&z->data->dmn,&x->data->dmn);
		z->data->dmn.num += y->data->dmn.num;
		break;
	}
	return z;
}

LOCAL NODEP subtract(x,y)
NODEP x,y;
{
	NODEP z;

	switch (x->type) {
	case INODE:
		z = allocinode(x->data->lng - y->data->lng);
		break;
	case FNODE:
		z = allocfnode(x->data->frt.numer * y->data->frt.denom -
						y->data->frt.numer * x->data->frt.denom,
						x->data->frt.denom * y->data->frt.denom);
		break;
	case DNODE:
		z = allocdnode(x->data->dbl - y->data->dbl);
		break;
	case NNODE:
		z = allocnnode();
		if (!check_equ_dim(&x->data->dmn,&y->data->dmn))
			printf("Subtracting unlike dimensions!\n");
		else
			copydnum(&z->data->dmn,&x->data->dmn);
		z->data->dmn.num -= y->data->dmn.num;
		break;
	}
	return z;
}

LOCAL NODEP multiply(x,y)
NODEP x,y;
{
	NODEP z;
	int i;

	switch (x->type) {
	case INODE:
		z = allocinode(x->data->lng * y->data->lng);
		break;
	case FNODE:
		z = allocfnode(x->data->frt.numer * y->data->frt.numer,
						x->data->frt.denom * y->data->frt.denom);
		break;
	case DNODE:
		z = allocdnode(x->data->dbl * y->data->dbl);
		break;
	case NNODE:
		z = allocnnode();
		z->data->dmn.num = x->data->dmn.num * y->data->dmn.num;
		for (i=0; i<MAXDIM; i++)
			z->data->dmn.di[i] = x->data->dmn.di[i] + y->data->dmn.di[i];
		break;
	}
	return z;
}

LOCAL NODEP divide(x,y)	/* Should not recieve a list or string node */
NODEP x,y;			/* the node types of x & y must be the same */
{	/* This is the longest of the operator functions because it has to
		cope with possible division by zero */
	NODEP z;
	int i;

	switch (x->type) {
	case INODE:		/* division of two integers results in a double */
		if (y->data->lng==0L) {
			error = MTHERR;
			return NULL;
		}
		z = allocdnode((double) x->data->lng / (double) y->data->lng );
		break;
	case FNODE:
		if (y->data->frt.numer == 0) {
			error = MTHERR;
			return NULL;
		}
		z = allocfnode(x->data->frt.numer*y->data->frt.denom,
						y->data->frt.numer*x->data->frt.denom);
		break;
	case DNODE:
		if (y->data->dbl==0.0) {
			error = MTHERR;
			return NULL;
		}
		z = allocdnode(x->data->dbl / y->data->dbl);
		break;
	case NNODE:
		if (y->data->dmn.num==0.0) {
			error = MTHERR;
			return NULL;
		}
		z = allocnnode();
		z->data->dmn.num = x->data->dmn.num / y->data->dmn.num;
		for (i=0; i<MAXDIM; i++)
			z->data->dmn.di[i] = x->data->dmn.di[i] - y->data->dmn.di[i];
		break;
	}
	return z;	/* Calling procedure should deallocate x & y */
}

LOCAL NODEP power(x,y)	/* y may not be dimensioned */
NODEP x,y;
{
	NODEP z;
	int i;

	if (y->type==NNODE) {
		error = TYPERR;
		return NULL;
	}
	change_type(y,DNODE);		/* promote y to dbl */
	if (x->type!=NNODE) {
		change_type(x,DNODE);	/* promote x to dbl */
		z = allocdnode(pow(x->data->dbl,y->data->dbl));
	} else {
		z = allocnnode();
		z->data->dmn.num = pow(x->data->dmn.num,y->data->dbl);
		for (i=0; i<MAXDIM; i++)
			z->data->dmn.di[i] = (int) (x->data->dmn.di[i] * y->data->dbl);
	}
	return z;
}

LOCAL NODEP and(x,y)
NODEP x,y;
{
	if (x->type!=INODE || y->type!=INODE) {
		error = TYPERR;
		return NULL;
	}
	return allocinode(x->data->lng && y->data->lng);
}

LOCAL NODEP logical_and(x,y)
NODEP x,y;
{
	if (x->type!=INODE || y->type!=INODE) {
		error = TYPERR;
		return NULL;
	}
	return allocinode(x->data->lng & y->data->lng);
}

LOCAL NODEP not(x)
NODEP x;
{
	if (x->type != INODE) {
		error = TYPERR;
		return NULL;
	}
	return allocinode(!x->data->lng);
}

LOCAL NODEP logical_not(x)
NODEP x;
{
	if (x->type != INODE) {
		error = TYPERR;
		return NULL;
	}
	return allocinode(~x->data->lng);
}

LOCAL NODEP or(x,y)
NODEP x,y;
{
	if (x->type != INODE) {
		error = TYPERR;
		return NULL;
	}
	return allocinode(x->data->lng || y->data->lng);
}

LOCAL NODEP logical_or(x,y)
NODEP x,y;
{
	if (x->type != INODE) {
		error = TYPERR;
		return NULL;
	}
	return allocinode(x->data->lng | y->data->lng);
}

LOCAL NODEP modulo(x,y)
NODEP x,y;
{
	if (x->type != INODE || y->type != INODE) {
		error = TYPERR;
		return NULL;
	}
	return allocinode(x->data->lng % y->data->lng);
}

LOCAL NODEP minus(x)	/* Unary minus */
NODEP x;
{
	NODEP z;

	switch (x->type) {
	case INODE:
		return allocinode(- x->data->lng);
	case FNODE:
		return allocfnode(- x->data->frt.numer, x->data->frt.denom);
	case DNODE:
		return allocdnode(- x->data->dbl);
	case NNODE:
		z = allocnnode();
		copydnum(&z->data->dmn,&x->data->dmn);
		z->data->dmn.num = - z->data->dmn.num;
		return z;
	}
	return NULL;
}


/*-------- Operator data --------*/


EXPORT char const operator[OPS][OPLEN] = {
		"Null",	/* cannot match anything because of mixed upper/lower case */
		"+","-","*","/","^","AND","NOT","!","OR","MOD","-",
		"&&","&","~","||","|"
	};
EXPORT int const priority[OPS] = {
		0,
		120,120,130,130,135,50,140,132,40,130,140,
		50,80,140,40,60
	};
/* Operator type: 0=NULL, 1=infix (binary), 2=prefix, 3=postfix */
EXPORT int const optype[OPS] = {
		0,
		1,1,1,1,1,1,2,3,1,1,2,
		1,1,2,1,1
	};
EXPORT NODEP (*do_op[OPS])() = {	/* array of pointers to functions */
		add,								/* filler for 0 index */
		add,subtract,multiply,divide,power,and,not,factorial,or,modulo,minus,
		and, logical_and, logical_not, or, logical_or
	};