/* General purpose routines, initialisation & tables */
#include "evalint.h"
#include <string.h>
#include <math.h>
#include "support.h"

int eval_error;

/* List of functions with their names */
onefunc functions[] = {
    { sin,   "sin" },
    { cos,   "cos" },
    { tan,   "tan" },
    { asin,  "asin" },
    { acos,  "acos" },
    { atan,  "atan" },
    { sinh,  "sinh" },
    { cosh,  "cosh" },
    { tanh,  "tanh" },
    { asinh, "asinh" },
    { acosh, "acosh" },
    { atanh, "atanh" },
    { exp,   "exp" },
    { exp10, "exp10" },
    { fabs,  "abs" },
    { log,   "log" },
    { log10, "log10" },
    { sqrt,  "sqrt" },
    { sqr,   "sqr" },
    { gamma, "gamma" },
    { neg,   "-" },
};

/* Simplifications done by eval when requested (PAT flag) */
/* x and y can be any expression */
char *simplifications[] = {
    "x+(-y)",           "x-y",
    "-x+y",             "y-x",
    "x-(-y)",           "x+y",
    "-(-x)",            "x",
    "sqr(x^y)",         "x^(y*2)",
    "0+x",              "x",
    "x+0",              "x",
    "x-0",              "x",
    "0-x",              "-x",
    "x-x",              "0",
    "x*0",              "0",
    "0*x",              "0",
    "1*x",              "x",
    "x*1",              "x",
    "-1*x",             "-x",
    "x*-1",             "-x",
    "x/1",              "x",
    "-x/-y",            "x/y",
    "1^x",              "1",
    "x^0",              "1",
    "x^1",              "x",
    "x^2",              "sqr(x)",
    "x^-1",             "1/x",
    "sin(-x)",          "-sin(x)",
    "cos(-x)",          "cos(x)",
    "tan(atan(x))",     "x",
    "tan(-x)",          "-tan(x)",
    "sin(x)/cos(x)",    "tan(x)",
    "abs(abs(x))",      "abs(x)",
    "abs(-x)",          "abs(x)",
    "10^x",             "exp10(x)",
    "log(exp(x))",      "x",
    "log10(exp10(x))",  "x",
    "log10(exp(x))",    "x*log10(e)",
    "log(exp10(x))",    "x*log(10)",
    "sinh(asinh(x))",   "x",
    "asinh(sinh(x))",   "x",
/* More 'doubtful' simplifications: eg 0/y is 0 for y != 0, but is
  undefined for y = 0 */
    "0/x",              "0",
    "x/sqr(x)",         "1/x"
};

#define NBSIMP (sizeof(simplifications) / (2 * sizeof(char *)))

/* For storing the "compiled" version of the substitutions */
subst *simp[NBSIMP + 1];

/* Rules for derivation */
/* x and y can be any expression, dx is the differential of x, etc */
/* Don't use d<name> in the left hand column ... */
char *derivees[] = {
    "x+y",          "dx+dy",
    "x-y",          "dx-dy",
    "x*y",          "dx*y+x*dy",
    "x/y",          "(dx*y-x*dy)/sqr(y)",
    "x^y",          "y*x^(y-1)*dx+x^y*log(x)*dy",
    "-x",           "-dx",
    "sin(x)",       "cos(x)*dx",
    "cos(x)",       "-sin(x)*dx",
    "tan(x)",       "(1+sqr(tan(x)))*dx",
    "asin(x)",      "dx/sqrt(1-sqr(x))",
    "acos(x)",      "-dx/sqrt(1-sqr(x))",
    "atan(x)",      "dx/(1+sqr(x))",
    "sinh(x)",      "cosh(x)*dx",
    "cosh(x)",      "sinh(x)*dx",
    "tanh(x)",      "(1-sqr(tanh(x)))*dx",
    "asinh(x)",     "dx/sqrt(sqr(x)+1)",
    "acosh(x)",     "dx/sqrt(sqr(x)-1)",
    "atanh(x)",     "dx/(1-sqr(x))",
    "exp(x)",       "exp(x)*dx",
    "exp10(x)",     "exp10(x)*log(10)*dx",
    "log(x)",       "dx/x",
    "log10(x)",     "dx/(x*log(10))",
    "sqrt(x)",      "dx/(2*sqrt(x))",
    "sqr(x)",       "2*x*dx"
};

#define NBDERV (sizeof(derivees) / (2 * sizeof(char *)))

subst *derv[NBDERV + 1];

/* Initialise the evaluation system */
int init_expr()
{
    register int i;

    /* Prepare substitutions */
    for (i = 0; i < NBSIMP; i++)
	if (!(simp[i] = make_sub(simplifications[2 * i], simplifications[2 * i + 1])))
	    return(FALSE);

    simp[i] = NULL;

    /* prepare derivations */
    for (i = 0; i < NBDERV; i++)
	if (!(derv[i] = make_sub(derivees[2 * i], derivees[2 * i + 1])))
	    return(FALSE);

    derv[i] = NULL;

    return(TRUE);
}

/* clean up after use */
void cleanup_expr()
{
    register int i;

    for (i = 0; i < NBSIMP; i++)
	if (simp[i]) free_sub(simp[i]);

    for (i = 0; i < NBDERV; i++)
	if (derv[i]) free_sub(derv[i]);
}

/* Map /+-* errors to math.h exceptions */
double check_fp(res, res1, res2, op)
double res, res1, res2;
char *op;
{
    switch (_FPERR) {
	case FPEUND: return(except(UNDERFLOW, op, res1, res2, 0.0));
	case FPEOVF: return(except(OVERFLOW, op, res1, res2, sign(res1) * sign(res2) * HUGE));
	case FPEZDV: return(except(SING, op, res1, res2, sign(res1) * HUGE));
	case FPENAN: return(except(TLOSS, op, res1, res2, 0.0));
	case FPECOM: return(except(TLOSS, op, res1, res2, 0.0));
	default: return(res);
    }
}

/* Find functions address from name */
onefunc *locate_name(name)
char *name;
{
    register int i;
    register onefunc *f = NULL;

    for (i = 0; i < sizeof(functions) / sizeof(onefunc); i++)
	if (strcmp(functions[i].functext, name) == 0) { f = &functions[i]; break; }

    return(f);
}

/* Find function name from address */
onefunc *locate(func)
FUNCTION func;
{
    register int i;
    register onefunc *f = NULL;

    for (i = 0; i < sizeof(functions) / sizeof(onefunc); i++)
	if (functions[i].func == func) { f = &functions[i]; break; }

    return(f);
}