/**********
Copyright 1990 Regents of the University of California.  All rights reserved.
Author: 1987 Wayne A. Christopher, U. C. Berkeley CAD Group
**********/
//#define TRACE

#include "ngspice/ngspice.h"
#include "ngspice/ifsim.h"
#include "ngspice/iferrmsg.h"
#include "ngspice/inpdefs.h"
#include "ngspice/inpptree.h"
#include "inpxx.h"

#include "inpptree-parser.h"
#include "inpptree-parser-y.h"

#ifdef OLD_BISON
int PTparse(char **line, INPparseNode **p, CKTcircuit *ckt);
#endif

static INPparseNode *mkcon(double value);
static INPparseNode *mkb(int type, INPparseNode * left,
                         INPparseNode * right);
static INPparseNode *mkf(int type, INPparseNode * arg);
static int PTcheck(INPparseNode * p);
static INPparseNode *mkvnode(char *name);
static INPparseNode *mkinode(char *name);

static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum);

static void free_tree(INPparseNode *);
static void printTree(INPparseNode *);


/*
 * LAW for INPparseNode* generator and consumer functions:
 * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 *   Newly allocated structs shall be initialized with `usecnt' = 0
 *   When filling INPparseNode * slots of newly initialized structs
 *     their `usecnt' shall be incremented
 *   Generators pass the responsibility `to free' return values
 *     on to their invokers.
 *   Functions generally process args with exactly one of:
 *     - inc_usage(arg) if they insert an argument into a struct
 *     - release_tree(arg) if they don't make any use of it
 *     - pass it on to another function()
 *   Functions use the the result of a function invocations with one of:
 *     - inc_usage(result) if they insert the result into a struct
 *     - release_tree(result) if they don't make any use of it
 *     - pass it on to another function()
 *     - simply return the result
 *
 * mkfirst(first, second)
 *   is used to safely release its second argument,
 *   and return its first
 *
 */


static inline INPparseNode *
inc_usage(INPparseNode *p)
{
    if(p)
        p->usecnt ++;
    return p;
}


static void
dec_usage(INPparseNode *p)
{
    if(p && --p->usecnt <= 0)
        free_tree(p);
}


static void
release_tree(INPparseNode *p)
{
    if(p && p->usecnt <= 0)
        free_tree(p);
}


static INPparseNode *
mkfirst(INPparseNode *fst, INPparseNode *snd)
{
    if(fst) {
        fst->usecnt ++;
        release_tree(snd);
        fst->usecnt --;
    } else {
        release_tree(snd);
    }

    return fst;
}


static IFvalue *values = NULL;
static int *types;
static int numvalues;
static CKTcircuit *circuit;
static INPtables *tables;

extern IFsimulator *ft_sim;        /* XXX */

/* Some tables that the parser uses. */

static struct op {
    int number;
    char *name;
    void (*funcptr)(void);
} ops[] = {
    {
    PT_COMMA,  ",", NULL}, {
    PT_PLUS,   "+", (void(*)(void)) PTplus}, {
    PT_MINUS,  "-", (void(*)(void)) PTminus}, {
    PT_TIMES,  "*", (void(*)(void)) PTtimes}, {
    PT_DIVIDE, "/", (void(*)(void)) PTdivide}, {
    PT_POWER,  "^", (void(*)(void)) PTpower}
};

#define NUM_OPS (int)NUMELEMS(ops)

static struct func {
    char *name;
    int number;
    void (*funcptr)(void);
} funcs[] = {
    { "abs",    PTF_ABS,    (void(*)(void)) PTabs } ,
    { "acos",   PTF_ACOS,   (void(*)(void)) PTacos } ,
    { "acosh",  PTF_ACOSH,  (void(*)(void)) PTacosh } ,
    { "asin",   PTF_ASIN,   (void(*)(void)) PTasin } ,
    { "asinh",  PTF_ASINH,  (void(*)(void)) PTasinh } ,
    { "atan",   PTF_ATAN,   (void(*)(void)) PTatan } ,
    { "atanh",  PTF_ATANH,  (void(*)(void)) PTatanh } ,
    { "cos",    PTF_COS,    (void(*)(void)) PTcos } ,
    { "cosh",   PTF_COSH,   (void(*)(void)) PTcosh } ,
    { "exp",    PTF_EXP,    (void(*)(void)) PTexp } ,
    { "ln",     PTF_LOG,    (void(*)(void)) PTlog } ,
    { "log",    PTF_LOG,    (void(*)(void)) PTlog } ,
    { "log10",  PTF_LOG10,  (void(*)(void)) PTlog10 } ,
    { "sgn",    PTF_SGN,    (void(*)(void)) PTsgn } ,
    { "sin",    PTF_SIN,    (void(*)(void)) PTsin } ,
    { "sinh",   PTF_SINH,   (void(*)(void)) PTsinh } ,
    { "sqrt",   PTF_SQRT,   (void(*)(void)) PTsqrt } ,
    { "tan",    PTF_TAN,    (void(*)(void)) PTtan } ,
    { "tanh",   PTF_TANH,   (void(*)(void)) PTtanh } ,
    { "u",      PTF_USTEP,  (void(*)(void)) PTustep } ,
    { "uramp",  PTF_URAMP,  (void(*)(void)) PTuramp } ,
    { "ceil",   PTF_CEIL,   (void(*)(void)) PTceil } ,
    { "floor",  PTF_FLOOR,  (void(*)(void)) PTfloor } ,
    { "nint",   PTF_NINT,   (void(*)(void)) PTnint } ,
    { "-",      PTF_UMINUS, (void(*)(void)) PTuminus },
    /* MW. cif function added */
    { "u2",     PTF_USTEP2, (void(*)(void)) PTustep2},
    { "pwl",    PTF_PWL,    (void(*)(void)) PTpwl},
    { "pwl_derivative", PTF_PWL_DERIVATIVE, (void(*)(void)) PTpwl_derivative},
    { "eq0",    PTF_EQ0,    (void(*)(void)) PTeq0},
    { "ne0",    PTF_NE0,    (void(*)(void)) PTne0},
    { "gt0",    PTF_GT0,    (void(*)(void)) PTgt0},
    { "lt0",    PTF_LT0,    (void(*)(void)) PTlt0},
    { "ge0",    PTF_GE0,    (void(*)(void)) PTge0},
    { "le0",    PTF_LE0,    (void(*)(void)) PTle0},
    { "pow",    PTF_POW,    (void(*)(void)) PTpower},
    { "pwr",    PTF_PWR,    (void(*)(void)) PTpwr},
    { "min",    PTF_MIN,    (void(*)(void)) PTmin},
    { "max",    PTF_MAX,    (void(*)(void)) PTmax},
} ;

#define NUM_FUNCS (int)NUMELEMS(funcs)

/* These are all the constants any sane person needs. */

static struct constant {
    char *name;
    double value;
} constants[] = {
    {
    "e", M_E}, {
    "pi", M_PI}
};

#define NUM_CONSTANTS (int)NUMELEMS(constants)

/* Parse the expression in *line as far as possible, and return the parse
 * tree obtained.  If there is an error, *pt will be set to NULL and an error
 * message will be printed.
 */

void
INPgetTree(char **line, INPparseTree ** pt, CKTcircuit *ckt, INPtables * tab)
{
    INPparseNode *p = NULL;
    int i, rv;

    values = NULL;
    types = NULL;
    numvalues = 0;

    circuit = ckt;
    tables = tab;

#ifdef TRACE
    fprintf(stderr,"%s, line = \"%s\"\n", __func__, *line);
#endif

    rv = PTparse(line, &p, ckt);

    if (rv || !p || !PTcheck(p)) {

        *pt = NULL;
        release_tree(p);

    } else {

        (*pt) = TMALLOC(INPparseTree, 1);

        (*pt)->p.numVars = numvalues;
        (*pt)->p.varTypes = types;
        (*pt)->p.vars = values;
        (*pt)->p.IFeval = IFeval;
        (*pt)->tree = inc_usage(p);

        (*pt)->derivs = TMALLOC(INPparseNode *, numvalues);

        for (i = 0; i < numvalues; i++)
            (*pt)->derivs[i] = inc_usage(PTdifferentiate(p, i));

    }

    values = NULL;
    types = NULL;
    numvalues = 0;

    circuit = NULL;
    tables = NULL;
}

/* This routine takes the partial derivative of the parse tree with respect to
 * the i'th variable.  We try to do optimizations like getting rid of 0-valued
 * terms.
 *
 *** Note that in the interests of simplicity we share some subtrees between
 *** the function and its derivatives.  This means that you can't free the
 *** trees.
 */

static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum)
{
    INPparseNode *arg1 = NULL, *arg2 = NULL, *newp = NULL;

    switch (p->type) {
    case PT_TIME:
    case PT_TEMPERATURE:
    case PT_FREQUENCY:
    case PT_CONSTANT:
        newp = mkcon(0.0);
        break;

    case PT_VAR:
        /* Is this the variable we're differentiating wrt? */
        if (p->valueIndex == varnum)
            newp = mkcon(1.0);
        else
            newp = mkcon(0.0);
        break;

    case PT_PLUS:
    case PT_MINUS:
        arg1 = PTdifferentiate(p->left, varnum);
        arg2 = PTdifferentiate(p->right, varnum);
        newp = mkb(p->type, arg1, arg2);
        break;

    case PT_TIMES:
        /* d(a * b) = d(a) * b + d(b) * a */
        arg1 = PTdifferentiate(p->left, varnum);
        arg2 = PTdifferentiate(p->right, varnum);

        newp = mkb(PT_PLUS, mkb(PT_TIMES, arg1, p->right),
                   mkb(PT_TIMES, p->left, arg2));
        break;

    case PT_DIVIDE:
        /* d(a / b) = (d(a) * b - d(b) * a) / b^2 */
        arg1 = PTdifferentiate(p->left, varnum);
        arg2 = PTdifferentiate(p->right, varnum);

        newp = mkb(PT_DIVIDE, mkb(PT_MINUS, mkb(PT_TIMES, arg1,
                                                p->right), mkb(PT_TIMES,
                                                               p->left,
                                                               arg2)),
                   mkb(PT_POWER, p->right, mkcon(2.0)));
        break;

    case PT_POWER:
#define a  p->left
#define b  p->right
        if (b->type == PT_CONSTANT) {
            /*
             * D(a^C) = C * a^(C-1) * D(a)
             */
            arg1 = PTdifferentiate(a, varnum);

            newp = mkb(PT_TIMES, mkb(PT_TIMES,
                                     mkcon(b->constant),
                                     mkb(PT_POWER, a,
                                         mkcon(b->constant - 1))),
                       arg1);
        } else {
            /*
             * D(a^b) = D(exp(b*log(a)))
             *              = exp(b*log(a)) * D(b*log(a))
             *              = exp(b*log(a)) * (D(b)*log(a) + b*D(a)/a)
             */
            arg1 = PTdifferentiate(a, varnum);
            arg2 = PTdifferentiate(b, varnum);
            newp = mkb(PT_TIMES, mkf(PTF_EXP, mkb(PT_TIMES,
                                                  b, mkf(PTF_LOG,
                                                         a))),
                       mkb(PT_PLUS,
                           mkb(PT_TIMES, b,
                               mkb(PT_DIVIDE, arg1, a)),
                           mkb(PT_TIMES, arg2, mkf(PTF_LOG, a))));
        }
#undef b
#undef a
        break;

    case PT_TERN: /* ternary_fcn(cond,exp1,exp2) */
      // naive:
      //   d/d ternary_fcn(cond,exp1,exp2) --> ternary_fcn(cond, d/d exp1, d/d exp2)
      {
//        extern void printTree(INPparseNode *);
//
//        printf("debug: %s, PT_TERN: ", __func__);
//        printTree(p);
//        printf("\n");

        newp = mkb(PT_TERN, p->left, mkb(PT_COMMA,
                                         PTdifferentiate(p->right->left, varnum),
                                         PTdifferentiate(p->right->right, varnum)));

//        printf("debug, %s, returns; ", __func__);
//        printTree(newp);
//        printf("\n");

        return mkfirst(newp, p);
      }

    case PT_FUNCTION:
        /* Many cases.  Set arg1 to the derivative of the function,
         * and arg2 to the derivative of the argument.
         */
        switch (p->funcnum) {
        case PTF_ABS:                /* sgn(u) */
            arg1 = mkf(PTF_SGN, p->left);
            break;

        case PTF_SGN:
            arg1 = mkcon(0.0);
            break;

        case PTF_ACOS:                /* - 1 / sqrt(1 - u^2) */
            arg1 = mkb(PT_DIVIDE, mkcon(-1.0), mkf(PTF_SQRT,
                                                          mkb(PT_MINUS,
                                                              mkcon(1.0),
                                                              mkb(PT_POWER,
                                                                  p->left,
                                                                  mkcon(2.0)))));
            break;

        case PTF_ACOSH:        /* 1 / sqrt(u^2 - 1) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkf(PTF_SQRT,
                                                         mkb(PT_MINUS,
                                                             mkb(PT_POWER,
                                                                 p->left,
                                                                 mkcon(2.0)),
                                                             mkcon(1.0))));

            break;

        case PTF_ASIN:                /* 1 / sqrt(1 - u^2) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkf(PTF_SQRT,
                                                         mkb(PT_MINUS,
                                                             mkcon(1.0),
                                                             mkb(PT_POWER,
                                                                 p->left,
                                                                 mkcon(2.0)))));
            break;

        case PTF_ASINH:        /* 1 / sqrt(u^2 + 1) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkf(PTF_SQRT,
                                                         mkb(PT_PLUS,
                                                             mkb(PT_POWER,
                                                                 p->left,
                                                                 mkcon(2.0)),
                                                             mkcon(1.0))));
            break;

        case PTF_ATAN:                /* 1 / (1 + u^2) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkb(PT_PLUS,
                                                         mkb(PT_POWER,
                                                             p->left,
                                                             mkcon(2.0)),
                                                         mkcon(1.0)));
            break;

        case PTF_ATANH:        /* 1 / (1 - u^2) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkb(PT_MINUS,
                                                         mkcon(1.0),
                                                         mkb(PT_POWER,
                                                             p->left,
                                                             mkcon(2.0))));
            break;

        case PTF_COS:                /* - sin(u) */
            arg1 = mkf(PTF_UMINUS, mkf(PTF_SIN, p->left));
            break;

        case PTF_COSH:                /* sinh(u) */
            arg1 = mkf(PTF_SINH, p->left);
            break;

        case PTF_EXP:                /* exp(u) */
            arg1 = mkf(PTF_EXP, p->left);
            break;

        case PTF_LOG:               /* 1 / u */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), p->left);
            break;

        case PTF_LOG10:              /* log(e) / u */
            arg1 = mkb(PT_DIVIDE, mkcon(M_LOG10E), p->left);
            break;

        case PTF_SIN:                /* cos(u) */
            arg1 = mkf(PTF_COS, p->left);
            break;

        case PTF_SINH:                /* cosh(u) */
            arg1 = mkf(PTF_COSH, p->left);
            break;

        case PTF_SQRT:                /* 1 / (2 * sqrt(u)) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkb(PT_TIMES,
                                                         mkcon(2.0),
                                                         mkf(PTF_SQRT,
                                                             p->left)));
            break;

        case PTF_TAN:                /* 1 / (cos(u) ^ 2) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkb(PT_POWER,
                                                         mkf(PTF_COS,
                                                             p->left),
                                                         mkcon(2.0)));
            break;

        case PTF_TANH:                /* 1 / (cosh(u) ^ 2) */
            arg1 = mkb(PT_DIVIDE, mkcon(1.0), mkb(PT_POWER,
                                                         mkf(PTF_COSH,
                                                             p->left),
                                                         mkcon(2.0)));
            break;

        case PTF_USTEP:
        case PTF_EQ0:
        case PTF_NE0:
        case PTF_GT0:
        case PTF_LT0:
        case PTF_GE0:
        case PTF_LE0:
            arg1 = mkcon(0.0);
            break;

        case PTF_URAMP:
            arg1 = mkf(PTF_USTEP, p->left);
            break;

        case PTF_FLOOR:                /* naive: D(floor(u)) = 0 */
            arg1 = mkcon(0.0);
            break;

        case PTF_CEIL:                /* naive: D(ceil(u)) = 0 */
            arg1 = mkcon(0.0);
            break;

        case PTF_NINT:                /* naive: D(nint(u)) = 0 */
            arg1 = mkcon(0.0);
            break;


            /* MW. PTF_CIF for new cif function */
        case PTF_USTEP2: /* ustep2=uramp(x)-uramp(x-1) ustep2'=ustep(x)-ustep(x-1) */
            arg1 = mkb(PT_MINUS,
                       mkf(PTF_USTEP, p->left),
                       mkf(PTF_USTEP,
                           mkb(PT_MINUS,
                               p->left,
                               mkcon(1.0))));
            break;

        case PTF_UMINUS:    /* - 1 ; like a constant (was 0 !) */
            arg1 = mkcon(-1.0);
            break;

        case PTF_PWL: /* PWL(var, x1, y1, x2, y2, ... a const list) */
            arg1 = mkf(PTF_PWL_DERIVATIVE, p->left);
            arg1->data = p->data;
            break;

        case PTF_PWL_DERIVATIVE: /* d/dvar PWL(var, ...) */
            arg1 = mkcon(0.0);
            break;

        case PTF_MIN:
        case PTF_MAX:
        /* min(a,b) -->   (a<b)       ? a : b
        *           -->   ((a-b) < 0) ? a : b
        */
        {
            INPparseNode *a = p->left->left;
            INPparseNode *b = p->left->right;
            int comparison = (p->funcnum == PTF_MIN) ? PTF_LT0 : PTF_GT0;
#ifdef TRACE
            printf("debug: %s, PTF_MIN: ", __func__);
            printTree(p);
            printf("\n");
            printf("debug: %s, PTF_MIN, a: ", __func__);
            printTree(a);
            printf("\n");
            printf("debug: %s, PTF_MIN, b: ", __func__);
            printTree(b);
            printf("\n");
#endif
            newp = mkb(PT_TERN,
                       mkf(comparison, mkb(PT_MINUS, a, b)),
                       mkb(PT_COMMA,
                           PTdifferentiate(a, varnum),
                           PTdifferentiate(b, varnum)));
#ifdef TRACE
            printf("debug, %s, returns; ", __func__);
            printTree(newp);
            printf("\n");
#endif
            return mkfirst(newp, p);
        }

        break;

        case PTF_POW:
        {
        /*
        pow(a,b)
        p->left: ','    p->left->left: a       p->left->right: b
        */
#define a  p->left->left
#define b  p->left->right

            if (b->type == PT_CONSTANT) {
                /*
                 * D(f^C) = C * f^(C-1) * D(f)
                 */
                arg1 = PTdifferentiate(a, varnum);

                newp = mkb(PT_TIMES, mkb(PT_TIMES,
                                     mkcon(b->constant),
                                     mkb(PT_POWER, a,
                                         mkcon(b->constant - 1))),
                             arg1);
#ifdef TRACE
            printf("pow, %s, returns; ", __func__);
            printTree(newp);
            printf("\n");
#endif
            } else {
            /*
             * D(f^g) = D(exp(g*log(f)))
             *              = exp(g*log(f)) * D(g*log(f))
             *              = exp(g*log(f)) * (D(g)*log(f) + g*D(f)/f)
             */
             arg1 = PTdifferentiate(a, varnum);
             arg2 = PTdifferentiate(b, varnum);
             newp = mkb(PT_TIMES, mkf(PTF_EXP, mkb(PT_TIMES,
                                                b, mkf(PTF_LOG,
                                                a))),
                                mkb(PT_PLUS,
                                    mkb(PT_TIMES, b,
                                    mkb(PT_DIVIDE, arg1, a)),
                                    mkb(PT_TIMES, arg2, mkf(PTF_LOG, a))));
            }
            return mkfirst(newp, p);
#undef b
#undef a
        }

        break;

        case PTF_PWR:
        {
        /*
        pwr(a,b)
        p->left: ','    p->left->left: a       p->left->right: b
        */
#define a  p->left->left
#define b  p->left->right
            if (b->type == PT_CONSTANT) {
                /* b is a constant
                 *
                 * f(a,b) = signum(a) * abs(a)^b
                 *        = signum(a) * exp(b*log(abs(a)))
                 * D(f)   = signum(a) * D(exp(b*log(abs(a))))
                 *        = signum(a) * exp(b*log(abs(a))) * D(b*log(abs(a)))
                 *        = signum(a) * abs(a)^b * D(b*log(abs(a)))
                 *        = signum(a) * abs(a)^b * b * 1/abs(a) * D(abs(a))
                 *        = signum(a) * abs(a)^(b-1) * b * D(abs(a))
                 *        = signum(a) * abs(a)^(b-1) * b * signum(a) * D(a)
                 *        = abs(a)^(b-1) * b * D(a)
                 */
                arg1 = PTdifferentiate(a, varnum);

                newp = mkb(PT_TIMES,
                           mkb(PT_TIMES,
                               mkcon(b->constant),
                               mkb(PT_POWER,
                                   mkf(PTF_ABS, a),
                                   mkcon(b->constant - 1.0))),
                           arg1);
#ifdef TRACE
            printf("pwr, %s, returns; ", __func__);
            printTree(newp);
            printf("\n");
#endif
            } else {
            /* b is a function
             *
             * f(a,b) = signum(a) * abs(a)^b
             *        = signum(a) * exp(b*log(abs(a)))
             * D(f)   = signum(a) * D(exp(b*log(abs(a))))
             *        = signum(a) * exp(b*log(abs(a))) * D(b*log(abs(a)))
             *        = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b * D(log(abs(a))))
             *        = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b * 1/abs(a) * D(abs(a)))
             *        = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b * 1/abs(a) * signum(a)*D(a))
             *        = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b/a*D(a))
             *        = signum(a) * exp(b*log(abs(a))) * D(b) * log(abs(a) + signum(a) * exp(b*log(abs(a))) / a * b * D(a)
             *        = signum(a) * exp(b*log(abs(a))) * D(b) * log(abs(a) + abs(a)^(b-1) * b * D(a)
             */
             arg1 = PTdifferentiate(a, varnum);
             arg2 = PTdifferentiate(b, varnum);
             newp = mkb(PT_PLUS,
                        mkb(PT_TIMES,
                            mkf(PTF_SGN, a),
                            mkb(PT_TIMES,
                                mkb(PT_POWER, mkf(PTF_ABS, a), b),
                                mkb(PT_TIMES, arg2,
                                    mkf(PTF_LOG, mkf(PTF_ABS, a))))),
                        mkb(PT_TIMES,
                            mkb(PT_TIMES,
                                mkb(PT_POWER,
                                    mkf(PTF_ABS, a),
                                    mkb(PT_MINUS, b, mkcon(1.0))),
                                b),
                            arg1));
            }
            return mkfirst(newp, p);
#undef b
#undef a
        }

        default:
            fprintf(stderr, "Internal Error: bad function # %d\n",
                    p->funcnum);
            return mkfirst(NULL, p);
        }

        arg2 = PTdifferentiate(p->left, varnum);

        newp = mkb(PT_TIMES, arg1, arg2);

        break;

    default:
        fprintf(stderr, "Internal error: bad node type %d\n", p->type);
        newp = NULL;
        break;
    }

    return mkfirst(newp, p);
}

static INPparseNode *mkcon(double value)
{
    INPparseNode *p = TMALLOC(INPparseNode, 1);

    p->type = PT_CONSTANT;
    p->constant = value;
    p->usecnt = 0;

    return (p);
}

static INPparseNode *mkb(int type, INPparseNode * left,
                         INPparseNode * right)
{
    INPparseNode *p;
    int i;

    if ((right->type == PT_CONSTANT) && (left->type == PT_CONSTANT)) {
        double value;
        switch (type) {
        case PT_TIMES:
            value = left->constant * right->constant;
            return mkfirst(mkcon(value), mkfirst(left, right));

        case PT_DIVIDE:
            value = left->constant / right->constant;
            return mkfirst(mkcon(value), mkfirst(left, right));

        case PT_PLUS:
            value = left->constant + right->constant;
            return mkfirst(mkcon(value), mkfirst(left, right));

        case PT_MINUS:
            value = left->constant - right->constant;
            return mkfirst(mkcon(value), mkfirst(left, right));

        case PT_POWER:
            value = pow(left->constant, right->constant);
            return mkfirst(mkcon(value), mkfirst(left, right));
        }
    }
    switch (type) {
    case PT_TIMES:
        if ((left->type == PT_CONSTANT) && (left->constant == 0))
            return mkfirst(left, right);
        else if ((right->type == PT_CONSTANT) && (right->constant == 0))
            return mkfirst(right, left);
        else if ((left->type == PT_CONSTANT) && (left->constant == 1))
            return mkfirst(right, left);
        else if ((right->type == PT_CONSTANT) && (right->constant == 1))
            return mkfirst(left, right);
        break;

    case PT_DIVIDE:
        if ((left->type == PT_CONSTANT) && (left->constant == 0))
            return mkfirst(left, right);
        else if ((right->type == PT_CONSTANT) && (right->constant == 1))
            return mkfirst(left, right);
        break;

    case PT_PLUS:
        if ((left->type == PT_CONSTANT) && (left->constant == 0))
            return mkfirst(right, left);
        else if ((right->type == PT_CONSTANT) && (right->constant == 0))
            return mkfirst(left, right);
        break;

    case PT_MINUS:
        if ((right->type == PT_CONSTANT) && (right->constant == 0))
            return mkfirst(left, right);
        else if ((left->type == PT_CONSTANT) && (left->constant == 0))
            return mkfirst(mkf(PTF_UMINUS, right), left);
        break;

    case PT_POWER:
        if (right->type == PT_CONSTANT) {
            if (right->constant == 0)
                return mkfirst(mkcon(1.0), mkfirst(left, right));
            else if (right->constant == 1)
                return mkfirst(left, right);
        }
        break;

    case PT_TERN:
        if (left->type == PT_CONSTANT) {
            /*FIXME > 0.0, >= 0.5, != 0.0 or what ? */
            p = (left->constant != 0.0) ? right->left : right->right;
            return mkfirst(p, mkfirst(right, left));
        }
        if((right->left->type == PT_CONSTANT) &&
           (right->right->type == PT_CONSTANT) &&
           (right->left->constant == right->right->constant))
            return mkfirst(right->left, mkfirst(right, left));
        break;
     }

     p = TMALLOC(INPparseNode, 1);

     p->type = type;
     p->usecnt = 0;

     p->left = inc_usage(left);
     p->right = inc_usage(right);

    if(type == PT_TERN) {
        p->function = NULL;
        p->funcname = NULL;
        return (p);
    }


    for (i = 0; i < NUM_OPS; i++)
        if (ops[i].number == type)
            break;
    if (i == NUM_OPS) {
        fprintf(stderr, "Internal Error: bad type %d\n", type);
        return (NULL);
    }
    p->function = ops[i].funcptr;
    p->funcname = ops[i].name;

    return (p);
}

static INPparseNode *mkf(int type, INPparseNode * arg)
{
    INPparseNode *p;
    int i;

    for (i = 0; i < NUM_FUNCS; i++)
        if (funcs[i].number == type)
            break;
    if (i == NUM_FUNCS) {
        fprintf(stderr, "Internal Error: bad type %d\n", type);
        return (NULL);
    }

    if (arg->type == PT_CONSTANT) {
        double constval = PTunary(funcs[i].funcptr) (arg->constant);
        return mkfirst(mkcon(constval), arg);
    }

    p = TMALLOC(INPparseNode, 1);

    p->type = PT_FUNCTION;
    p->usecnt = 0;

    p->left = inc_usage(arg);

    p->funcnum = funcs[i].number;
    p->function = funcs[i].funcptr;
    p->funcname = funcs[i].name;

    p->data = NULL;

    return (p);
}

/* Check for remaining PT_PLACEHOLDERs in the parse tree.  Returns 1 if ok. */

static int PTcheck(INPparseNode * p)
{
    switch (p->type) {
    case PT_PLACEHOLDER:
        return (0);

    case PT_TIME:
    case PT_TEMPERATURE:
    case PT_FREQUENCY:
    case PT_CONSTANT:
    case PT_VAR:
        return (1);

    case PT_FUNCTION:
        return (PTcheck(p->left));

    case PT_PLUS:
    case PT_MINUS:
    case PT_TIMES:
    case PT_DIVIDE:
    case PT_POWER:
    case PT_COMMA:
        return (PTcheck(p->left) && PTcheck(p->right));
    case PT_TERN:
        return (PTcheck(p->left) && PTcheck(p->right->left) && PTcheck(p->right->right));

    default:
        fprintf(stderr, "Internal error: bad node type %d\n", p->type);
        return (0);
    }
}

/* Binop node. */

INPparseNode *PT_mkbnode(const char *opstr, INPparseNode * arg1,
                             INPparseNode * arg2)
{
    INPparseNode *p;
    int i;

    for (i = 0; i < NUM_OPS; i++)
        if (!strcmp(ops[i].name, opstr))
            break;

    if (i == NUM_OPS) {
        fprintf(stderr, "Internal Error: no such op num %s\n", opstr);
        return mkfirst(NULL, mkfirst(arg1, arg2));
    }

    p = TMALLOC(INPparseNode, 1);

    p->type = ops[i].number;
    p->usecnt = 0;

    p->funcname = ops[i].name;
    p->function = ops[i].funcptr;
    p->left = inc_usage(arg1);
    p->right = inc_usage(arg2);

    return (p);
}

/*
 * prepare_PTF_PWL()
 *   for the PWL(expr, points...) function
 *     collect the given points, which are expected to be given
 *       literal constant
 *   strip them from the INPparseNode
 *     and pass them as an opaque struct alongside the
 *     INPparseNode for the PWL(expr) function call
 *
 * Note:
 *  the functionINPgetTree() is missing a recursive decending simplifier
 *  as a consequence we can meet a PTF_UMINUS->PTF_CONSTANT
 *    instead of a plain PTF_CONSTANT here
 *  of course this can get arbitrarily more complex
 *    for example PTF_TIMES -> PTF_CONSTANT, PTF_CONSTANT
 *      etc.
 *  currently we support only PFT_CONST and PTF_UMINUS->PTF_CONST
 */

#define Breakpoint do { __asm__ __volatile__ ("int $03"); } while(0)

static INPparseNode *prepare_PTF_PWL(INPparseNode *p)
{
    INPparseNode *w;
    struct pwldata { int n; double *vals; } *data;
    int i;

    if (p->funcnum != PTF_PWL) {
        fprintf(stderr, "PWL-INFO: %s, very unexpected\n", __func__);
        controlled_exit(1);
    }

#ifdef TRACE
    fprintf(stderr, "PWL-INFO: %s  building a PTF_PWL\n", __func__);
#endif
    i = 0;
    for(w = p->left; w->type == PT_COMMA; w = w->left)
        i++;

    if (i<2 || (i%1)) {
        fprintf(stderr, "Error: PWL(expr, points...) needs an even and >=2 number of constant args\n");
        return mkfirst(NULL, p);
    }

    data = TMALLOC(struct pwldata, 1);
    data->vals = TMALLOC(double, i);

    data->n = i;

    p->data = (void *) data;

    for (w = p->left ; --i >= 0 ; w = w->left)
        if (w->right->type == PT_CONSTANT) {
            data->vals[i] = w->right->constant;
        } else if (w->right->type == PT_FUNCTION &&
                   w->right->funcnum == PTF_UMINUS &&
                   w->right->left->type == PT_CONSTANT) {
            data->vals[i] = - w->right->left->constant;
        } else {
            fprintf(stderr, "PWL-ERROR: %s, not a constant\n", __func__);
            fprintf(stderr, "   type = %d\n", w->right->type);
            //Breakpoint;
            fprintf(stderr, "Error: PWL(expr, points...) only *literal* points are supported\n");
            return mkfirst(NULL, p);
        }

#ifdef TRACE
    for (i = 0 ; i < data->n ; i += 2)
        fprintf(stderr, "  (%lf %lf)\n", data->vals[i], data->vals[i+1]);
#endif

    for (i = 2 ; i < data->n ; i += 2)
        if(data->vals[i-2] >= data->vals[i]) {
            fprintf(stderr, "Error: PWL(expr, points...) the abscissa of points must be ascending\n");
            return mkfirst(NULL, p);
        }

    /* strip all but the first arg,
     *   and attach the rest as opaque data to the INPparseNode
     */

    w = inc_usage(w);
    dec_usage(p->left);
    p->left = w;

    return (p);
}


INPparseNode *PT_mkfnode(const char *fname, INPparseNode * arg)
{
    int i;
    INPparseNode *p;
    char buf[128];

    /* Make sure the case is ok. */
    (void) strcpy(buf, fname);
    strtolower(buf);

    if(!strcmp("ternary_fcn", buf)) {

        if(arg->type == PT_COMMA && arg->left->type == PT_COMMA) {

            INPparseNode *arg1 = arg->left->left;
            INPparseNode *arg2 = arg->left->right;
            INPparseNode *arg3 = arg->right;

            p = TMALLOC(INPparseNode, 1);

            p->type = PT_TERN;
            p->usecnt = 0;

            p->left = inc_usage(arg1);
            p->right = inc_usage(mkb(PT_COMMA, arg2, arg3));

            return mkfirst(p, arg);
        }

        fprintf(stderr, "Error: bogus ternary_fcn form\n");
        return mkfirst(NULL, arg);
    }

    for (i = 0; i < NUM_FUNCS; i++)
        if (!strcmp(funcs[i].name, buf))
            break;

    if (i == NUM_FUNCS) {
        fprintf(stderr, "Error: no such function '%s'\n", buf);
        return mkfirst(NULL, arg);
    }

    p = TMALLOC(INPparseNode, 1);

    p->type = PT_FUNCTION;
    p->usecnt = 0;

    p->left = inc_usage(arg);
    p->funcname = funcs[i].name;
    p->funcnum = funcs[i].number;
    p->function = funcs[i].funcptr;
    p->data = NULL;

    if(p->funcnum == PTF_PWL)
        p = prepare_PTF_PWL(p);

    return (p);
}

static INPparseNode *mkvnode(char *name)
{
    INPparseNode *p = TMALLOC(INPparseNode, 1);

    int i;
    CKTnode *temp;

    INPtermInsert(circuit, &name, tables, &temp);
    for (i = 0; i < numvalues; i++)
        if ((types[i] == IF_NODE) && (values[i].nValue == temp))
            break;
    if (i == numvalues) {
        if (numvalues) {
            values = TREALLOC(IFvalue, values, numvalues + 1);
            types = TREALLOC(int, types, numvalues + 1);
        } else {
            values = TMALLOC(IFvalue, 1);
            types = TMALLOC(int, 1);
        }
        values[i].nValue = temp;
        types[i] = IF_NODE;
        numvalues++;
    }
    p->valueIndex = i;
    p->type = PT_VAR;
    p->usecnt = 0;

    return (p);
}

static INPparseNode *mkinode(char *name)
{
    INPparseNode *p = TMALLOC(INPparseNode, 1);

    int i;

    INPinsert(&name, tables);
    for (i = 0; i < numvalues; i++)
        if ((types[i] == IF_INSTANCE) && (values[i].uValue == name))
            break;
    if (i == numvalues) {
        if (numvalues) {
            values = TREALLOC(IFvalue, values, numvalues + 1);
            types = TREALLOC(int, types, numvalues + 1);
        } else {
            values = TMALLOC(IFvalue, 1);
            types = TMALLOC(int, 1);
        }
        values[i].uValue = name;
        types[i] = IF_INSTANCE;
        numvalues++;
    }
    p->valueIndex = i;
    p->type = PT_VAR;
    p->usecnt = 0;

    return (p);
}

/* Number node. */

INPparseNode *PT_mknnode(double number)
{
    struct INPparseNode *p;

    p = TMALLOC(INPparseNode, 1);

    p->type = PT_CONSTANT;
    p->usecnt = 0;

    p->constant = number;

    return (p);
}

/* String node. */

INPparseNode *PT_mksnode(const char *string, void *ckt)
{
    int i, j;
    char buf[128];
    INPparseNode *p;

    /* Make sure the case is ok. */
    (void) strcpy(buf, string);
    strtolower(buf);

    p = TMALLOC(INPparseNode, 1);

    p->usecnt = 0;

    if(!strcmp("time", buf)) {
        p->type = PT_TIME;
        p->data = ckt;
        return p;
    }

    if(!strcmp("temper", buf)) {
        p->type = PT_TEMPERATURE;
        p->data = ckt;
        return p;
    }

    if(!strcmp("hertz", buf)) {
        p->type = PT_FREQUENCY;
        p->data = ckt;
        return p;
    }

    /* First see if it's something special. */
    for (i = 0; i < ft_sim->numSpecSigs; i++)
        if (!strcmp(ft_sim->specSigs[i], buf))
            break;
    if (i < ft_sim->numSpecSigs) {
        for (j = 0; j < numvalues; j++)
            if ((types[j] == IF_STRING) && !strcmp(buf, values[i].sValue))
                break;
        if (j == numvalues) {
            if (numvalues) {
                values = TREALLOC(IFvalue, values, numvalues + 1);
                types = TREALLOC(int, types, numvalues + 1);
            } else {
                values = TMALLOC(IFvalue, 1);
                types = TMALLOC(int, 1);
            }
            values[i].sValue = TMALLOC(char, strlen(buf) + 1);
            strcpy(values[i].sValue, buf);
            types[i] = IF_STRING;
            numvalues++;
        }
        p->valueIndex = i;
        p->type = PT_VAR;
        return (p);
    }

    for (i = 0; i < NUM_CONSTANTS; i++)
        if (!strcmp(constants[i].name, buf))
            break;

    if (i == NUM_CONSTANTS) {
        /* We'd better save this in case it's part of i(something). */
        p->type = PT_PLACEHOLDER;
        p->funcname = copy(string);
    } else {
        p->type = PT_CONSTANT;
        p->constant = constants[i].value;
    }

    return (p);
}

/* The lexical analysis routine. */

int PTlex (YYSTYPE *lvalp, struct PTltype *llocp, char **line)
{
    double td;
    int err;
    static char *specials = " \t()^+-*/,";
    char *sbuf, *s;
    int token;

    sbuf = *line;
#ifdef TRACE
//    printf("entering lexer, sbuf = '%s', lastoken = %d, lasttype = %d\n",
//        sbuf, lasttoken, lasttype);
#endif
    while ((*sbuf == ' ') || (*sbuf == '\t'))
        sbuf++;

    llocp->start = sbuf;

    switch (*sbuf) {
    case '\0':
        token = 0;
        break;

    case '?':
    case ':':
    case ',':
    case '-':
    case '+':
    case '/':
    case '^':
    case '(':
    case ')':
        token = *sbuf++;
        break;

    case '*':
      if(sbuf[1] == '*') {
        sbuf += 2;
        token = '^';            /* `**' is exponentiation */
        break;
      } else {
        token = *sbuf++;
        break;
      }

    case '&':
      if(sbuf[1] == '&') {
        sbuf += 2;
        token = TOK_AND;
        break;
      } else {
        token = *sbuf++;
        break;
      }

    case '|':
      if(sbuf[1] == '|') {
        sbuf += 2;
        token = TOK_OR;
        break;
      } else {
        token = *sbuf++;
        break;
      }

    case '=':
      if(sbuf[1] == '=') {
        sbuf += 2;
        token = TOK_EQ;
        break;
      } else {
        token = *sbuf++;
        break;
      }

    case '!':
      if(sbuf[1] == '=') {
        sbuf += 2;
        token = TOK_NE;
        break;
      } else {
        token = *sbuf++;
        break;
      }

    case '>':
      if(sbuf[1] == '=') {
        sbuf += 2;
        token = TOK_GE;
        break;
      } else {
        sbuf += 1;
        token = TOK_GT;
        break;
      }

    case '<':
      if(sbuf[1] == '>') {
        sbuf += 2;
        token = TOK_NE;
        break;
      }
      else if(sbuf[1] == '=') {
        sbuf += 2;
        token = TOK_LE;
        break;
      } else {
        sbuf += 1;
        token = TOK_LT;
        break;
      }

    default:
        {
            int n1 = -1;
            int n2 = -1;
            int n3 = -1;
            int n4 = -1;
            int n  = -1;

            sscanf(sbuf, "%*1[vV] ( %n%*[^ \t,()]%n , %n%*[^ \t,()]%n )%n",
                   &n1, &n2, &n3, &n4, &n);
            if(n != -1) {
                token = TOK_pnode;
                lvalp->pnode = mkb(PT_MINUS,
                                   mkvnode(copy_substring(sbuf+n1, sbuf+n2)),
                                   mkvnode(copy_substring(sbuf+n3, sbuf+n4)));
                sbuf += n;
                break;
            }
        }

        {
            int n1 = -1;
            int n2 = -1;
            int n  = -1;

            sscanf(sbuf, "%*1[vV] ( %n%*[^ \t,()]%n )%n", &n1, &n2, &n);
            if(n != -1) {
                token = TOK_pnode;
                lvalp->pnode = mkvnode(copy_substring(sbuf+n1, sbuf+n2));
                sbuf += n;
                break;
            }
        }

        {
            int n1 = -1;
            int n2 = -1;
            int n  = -1;

            sscanf(sbuf, "%*1[iI] ( %n%*[^ \t,()]%n )%n", &n1, &n2, &n);
            if(n != -1) {
                token = TOK_pnode;
                lvalp->pnode = mkinode(copy_substring(sbuf+n1, sbuf+n2));
                sbuf += n;
                break;
            }
        }

        td = INPevaluate(&sbuf, &err, 1);
        if (err == OK) {
            token = TOK_NUM;
            lvalp->num = td;
        } else {
        char *tmp;
            token = TOK_STR;
            for (s = sbuf; *s; s++)
                if (strchr(specials, *s))
                    break;
            tmp = TMALLOC(char, s - sbuf + 1);
            strncpy(tmp, sbuf, (size_t) (s - sbuf));
            tmp[s - sbuf] = '\0';
            lvalp->str = tmp;
            sbuf = s;
        }
    }

    *line = sbuf;

#ifdef TRACE
//    printf("PTlexer: token = %d, type = %d, left = '%s'\n",
//        el.token, el.type, sbuf); */
#endif
    llocp->stop = sbuf;
    return (token);
}


void INPfreeTree(IFparseTree *ptree)
{
    INPparseTree *pt = (INPparseTree *) ptree;

    int i;

    if (!pt)
        return;

    for (i = 0; i < pt->p.numVars; i++)
        dec_usage(pt->derivs[i]);

    dec_usage(pt->tree);

    txfree(pt->derivs);
    txfree(pt->p.varTypes);
    txfree(pt->p.vars);
    txfree(pt);
}


void free_tree(INPparseNode *pt)
{
    if(!pt)
        return;

    if(pt->usecnt) {
        fprintf(stderr, "ERROR: fatal internal error, %s\n", __func__);
        controlled_exit(1);
    }

    switch (pt->type) {
    case PT_TIME:
    case PT_TEMPERATURE:
    case PT_FREQUENCY:
    case PT_CONSTANT:
    case PT_VAR:
        break;

    case PT_PLUS:
    case PT_MINUS:
    case PT_TIMES:
    case PT_DIVIDE:
    case PT_POWER:
    case PT_COMMA:
    case PT_TERN:
        dec_usage(pt->right);
    case PT_FUNCTION:
        dec_usage(pt->left);
        break;

    default:
        printf("oops");
        break;
    }

    if(pt->type == PT_FUNCTION && pt->funcnum == PTF_PWL) {
        struct pwldata { int n; double *vals; } *data = (struct pwldata*)(pt->data);
        if(data) {
            txfree(data->vals);
            txfree(data);
        }
    }

    txfree(pt);
}


/* Debugging stuff. */

void INPptPrint(char *str, IFparseTree * ptree)
{
    int i;

    printf("%s\n\t", str);
    printTree(((INPparseTree *) ptree)->tree);
    printf("\n");
    for (i = 0; i < ptree->numVars; i++) {
        printf("d / d v%d : ", i);
        printTree(((INPparseTree *) ptree)->derivs[i]);
        printf("\n");
    }
    return;
}

void printTree(INPparseNode * pt)
{
    switch (pt->type) {
    case PT_TIME:
        printf("time(ckt = %p)", pt->data);
        break;

    case PT_TEMPERATURE:
        printf("temperature(ckt = %p)", pt->data);
        break;

    case PT_FREQUENCY:
        printf("frequency(ckt = %p)", pt->data);
        break;

    case PT_CONSTANT:
        printf("%g", pt->constant);
        break;

    case PT_VAR:
        printf("v%d", pt->valueIndex);
        break;

    case PT_PLUS:
        printf("(");
        printTree(pt->left);
        printf(") + (");
        printTree(pt->right);
        printf(")");
        break;

    case PT_MINUS:
        printf("(");
        printTree(pt->left);
        printf(") - (");
        printTree(pt->right);
        printf(")");
        break;

    case PT_TIMES:
        printf("(");
        printTree(pt->left);
        printf(") * (");
        printTree(pt->right);
        printf(")");
        break;

    case PT_DIVIDE:
        printf("(");
        printTree(pt->left);
        printf(") / (");
        printTree(pt->right);
        printf(")");
        break;

    case PT_POWER:
        printf("(");
        printTree(pt->left);
        printf(") ^ (");
        printTree(pt->right);
        printf(")");
        break;

    case PT_COMMA:
        printf("(");
        printTree(pt->left);
        printf(") , (");
        printTree(pt->right);
        printf(")");
        break;

    case PT_FUNCTION:
        printf("%s (", pt->funcname);
        printTree(pt->left);
        printf(")");
        break;

    case PT_TERN:
        printf("ternary_fcn (");
        printTree(pt->left);
        printf(") , (");
        printTree(pt->right);
        printf(")");
        break;

    default:
        printf("oops");
        break;
    }
    return;
}

