inpptree.c

ngspice又一个电子CAD仿真软件代码.功能更全

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

#include "ngspice.h"
#include <stdio.h>
#include <ctype.h>
#include "ifsim.h"
#include "iferrmsg.h"
#include "inpdefs.h"
#include "inpptree.h"
#include "inp.h"

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 *PTparse(char **line);
static INPparseNode *makepnode(PTelement * elem);
static INPparseNode *mkbnode(int opnum, INPparseNode * arg1,
			     INPparseNode * arg2);
static INPparseNode *mkfnode(char *fname, INPparseNode * arg);
static INPparseNode *mknnode(double number);
static INPparseNode *mksnode(char *string);
static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum);
static PTelement *PTlexer(char **line);

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



extern IFsimulator *ft_sim;	/* XXX */

/* Some tables that the parser uses. */

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

#define NUM_OPS (sizeof (ops) / sizeof (struct op))

static struct func {
    char *name;
    int number;
    double (*funcptr) ();
} funcs[] = {
    { "abs",    PTF_ABS,    PTabs } ,
    { "acos",   PTF_ACOS,   PTacos } ,
    { "acosh",  PTF_ACOSH,  PTacosh } ,
    { "asin",   PTF_ASIN,   PTasin } ,
    { "asinh",  PTF_ASINH,  PTasinh } ,
    { "atan",   PTF_ATAN,   PTatan } ,
    { "atanh",  PTF_ATANH,  PTatanh } ,
    { "cos",    PTF_COS,    PTcos } ,
    { "cosh",   PTF_COSH,   PTcosh } ,
    { "exp",    PTF_EXP,    PTexp } ,
    { "ln",     PTF_LN,     PTln } ,
    { "log",    PTF_LOG,    PTlog } ,
    { "sgn",    PTF_SGN,    PTsgn } ,
    { "sin",    PTF_SIN,    PTsin } ,
    { "sinh",   PTF_SINH,   PTsinh } ,
    { "sqrt",   PTF_SQRT,   PTsqrt } ,
    { "tan",    PTF_TAN,    PTtan } ,
    { "tanh",   PTF_TANH,   PTtanh } ,
    { "u",   	PTF_USTEP,  PTustep } ,
    { "uramp",  PTF_URAMP,  PTuramp } ,
    { "-",      PTF_UMINUS, PTuminus },
    /* MW. cif function added */
    { "u2",	PTF_USTEP2, PTustep2}
} ;

#define NUM_FUNCS (sizeof (funcs) / sizeof (struct func))

/* 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 (sizeof (constants) / sizeof (struct constant))

/* 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, void *ckt, INPtables * tab)
{
    INPparseNode *p;
    int i;

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

    circuit = ckt;
    tables = tab;

    p = PTparse(line);

    if (!p || !PTcheck(p)) {
	*pt = NULL;
	return;
    }

    (*pt) = (INPparseTree *) MALLOC(sizeof(INPparseTree));

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

    (*pt)->derivs = (INPparseNode **)
	MALLOC(numvalues * sizeof(INPparseNode *));

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

    return;
}

/* 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, *newp;

/* printf("differentiating: "); printTree(p); printf(" wrt var %d\n", varnum);*/

    switch (p->type) {
    case PT_CONSTANT:
	newp = mkcon((double) 0);
	break;

    case PT_VAR:
	/* Is this the variable we're differentiating wrt? */
	if (p->valueIndex == varnum)
	    newp = mkcon((double) 1);
	else
	    newp = mkcon((double) 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((double) 2)));
	break;

    case PT_POWER:
	/* Two cases... If the power is a constant then we're cool.
	 * Otherwise we have to be tricky.
	 */
	if (p->right->type == PT_CONSTANT) {
	    arg1 = PTdifferentiate(p->left, varnum);

	    newp = mkb(PT_TIMES, mkb(PT_TIMES,
				     mkcon(p->right->constant),
				     mkb(PT_POWER, p->left,
					 mkcon(p->right->constant - 1))),
		       arg1);
	} else {
	    /* This is complicated.  f(x) ^ g(x) ->
	     * exp(y(x) * ln(f(x)) ...
	     */
	    arg1 = PTdifferentiate(p->left, varnum);
	    arg2 = PTdifferentiate(p->right, varnum);
	    newp = mkb(PT_TIMES, mkf(PTF_EXP, mkb(PT_TIMES,
						  p->right, mkf(PTF_LN,
								p->left))),
		       mkb(PT_PLUS,
			   mkb(PT_TIMES, p->right,
			       mkb(PT_DIVIDE, arg1, p->left)),
			   mkb(PT_TIMES, arg2, mkf(PTF_LN, /*arg1*/p->left))));
		                                        /*changed by HT, '05/06/29*/

	}
	break;

    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, 0); */
	    arg1 = mkf(PTF_SGN, p->left);
	    break;

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

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

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

	    break;

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

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

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

	case PTF_ATANH:	/* 1 / (1 - u^2) */
	    arg1 = mkb(PT_DIVIDE, mkcon((double) 1), mkb(PT_MINUS,
							 mkcon((double) 1),
							 mkb(PT_POWER,
							     p->left,
							     mkcon((double)
								   2))));
	    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, 0); */
	    arg1 = mkf(PTF_EXP, p->left);
	    break;

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

	case PTF_LOG:		/* log(e) / u */
	    arg1 = mkb(PT_DIVIDE, mkcon((double) 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((double) 1), mkb(PT_TIMES,
							 mkcon((double) 2),
							 mkf(PTF_SQRT,
							     p->left)));
	    break;

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

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

	case PTF_USTEP:
	    arg1 = mkcon((double) 0.0);
	    break;

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

	    /* MW. PTF_CIF for new cif function */
	case PTF_USTEP2:
	    arg1 = mkcon((double) 0.0);
	    break;
	    
        case PTF_UMINUS:    /* - 1 ; like a constant (was 0 !) */
            arg1 = mkcon((double) - 1.0);
            break;

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

	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;
    }

/* printf("result is: "); printTree(newp); printf("\n"); */
    return (newp);
}

static INPparseNode *mkcon(double value)
{
    INPparseNode *p = (INPparseNode *) MALLOC(sizeof(INPparseNode));

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

    return (p);
}

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

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

	case PT_DIVIDE:
	    return (mkcon(left->constant / right->constant));

	case PT_PLUS:
	    return (mkcon(left->constant + right->constant));

	case PT_MINUS:
	    return (mkcon(left->constant - right->constant));

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

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

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

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

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

    p->type = type;
    p->left = left;
    p->right = right;

    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 = (INPparseNode *) MALLOC(sizeof(INPparseNode));
    int i;
    double constval;

    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) {
	constval = ((*funcs[i].funcptr) (arg->constant));
	return (mkcon(constval));
    }

    p->type = PT_FUNCTION;
    p->left = arg;

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

    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_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:
	return (PTcheck(p->left) && PTcheck(p->right));

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

/* The operator-precedence table for the parser. */

#define G 1			/* Greater than. */