ascend/compiler/relation.c

/*  ASCEND modelling environment
    Copyright (C) 2006, 2010 Carnegie Mellon University
    Copyright (C) 1990, 1993, 1994 Thomas Guthrie Epperly
    Copyright (C) 1993, 1994, 1995  Kirk Andre' Abbott
    Copyright (C) 1996 Benjamin Andrew Allan

    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, 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, see <http://www.gnu.org/licenses/>.
*//*
    @file
    Relation construction routines
*//*
    by Tom Epperly
    Created: 1/30/90
    Last in CVS $Revision: 1.32 $ $Date: 1998/03/17 22:09:24 $ $Author: ballan $
*/

#include <math.h>
#include <stdarg.h>
#include <ascend/general/platform.h>
#include <ascend/general/ascMalloc.h>
#include <ascend/general/panic.h>
#include <ascend/general/pool.h>
#include <ascend/general/list.h>
/* #include <general/pairlist.h> */
#include <ascend/general/stack.h>
#include <ascend/general/dstring.h>

#include "expr_types.h"
#include "name.h"
#include "nameio.h"
#include "instance_enum.h"
#include "bintoken.h"
#include "exprs.h"
#include "exprio.h"
#include "value_type.h"
#include "evaluate.h"
#include "forvars.h"
#include "find.h"
#include "findpath.h"
#include "sets.h"
#include "setinstval.h"
#include "instance_io.h"
#include "extcall.h"
#include "relation_util.h"
#include "rel_common.h"
#include "rel_blackbox.h"
#include "temp.h"
#include "atomvalue.h"
#include "mathinst.h"
#include "instquery.h"
#include "tmpnum.h"
#include "vlist.h"
#include "relation.h"

/*
 * internal form of RelationRelop for lval or rval use.
 */
#define RelRelop(r) ((r)->share->s.relop)

#define SUM 1
#define PROD 0
#ifndef abs
#define abs(a) ( ((a)>0) ? (a) : (-(a)) )
#endif

/*
 * Some global and exported variables.
 */
struct gl_list_t *g_relation_var_list = NULL;

int g_simplify_relations = 1;

int g_ExternalNodeStamps=0; /* incremented each time an new external
               * statement is seen */

/* fwd declaration */
static union RelationTermUnion
*CopyRelationSide(union RelationTermUnion *, unsigned long);

#ifdef THIS_IS_AN_UNUSED_FUNCTION
static
unsigned long ExprLength(register CONST struct Expr *start,
             register CONST struct Expr *stop)
{
  register unsigned long result=0;
  while(start!=stop){
    start = NextExpr(start);
    result++;
  }
  return result;
}
#endif /* THIS_IS_AN_UNUSED_FUNCTION */


static
void FigureOutError(struct value_t value,
            enum relation_errors *err,
            enum find_errors *ferr)
{
  assert(ValueKind(value)==error_value);
  *err = find_error;
  switch(ErrorValue(value)){
  case type_conflict:
  case dimension_conflict:
  case incorrect_name:
  case incorrect_such_that:
  case empty_choice:
  case empty_intersection:
  case temporary_variable_reused:
    *ferr = impossible_instance;
    break;
  case undefined_value:
    *ferr = undefined_instance;
    break;
  case name_unfound:
    *ferr = unmade_instance;
    break;
  default:
    ASC_PANIC("Unknown error type in FigureOutError.\n");
    break;
  }
}

/*-----------------------------------------------------------------------------
  CREATION AND MANAGEMENT OF RELATION TERMS

  It is cheaper to create relation terms in arrays the size of
  the union than individually because of operating system overhead.

  Lookout, the tokens have unionized: next they'll want a raise.
*/

/*
 * The define POOL_ALLOCTERM is for people who are pulling terms out
 * of a pool and promise to return them immediately.
 */

static pool_store_t g_term_pool = NULL;
/* A pool_store for 1 expression.
 * It is expected that objective functions will cause the
 * largest expressions.
 * Each time an expression is completed, it will be copied
 * into an array which can be created already knowing
 * its proper size. The array will be naturally in postfix.
 */

#define POOL_ALLOCTERM A_TERM(pool_get_element(g_term_pool))
/* get a token. Token is the size of the RelationTermUnion. */
#ifdef NDEBUG
#define PTINIT(x)
#else
#define PTINIT(x) TermUnionInit(x)
#endif
#define POOL_RESET pool_clear_store(g_term_pool)
/* reset the pool for next expression */

#ifndef NDEBUG
/*
 * this function zeros a termunion ptr contents. tu must not be NULL.
 */
static void TermUnionInit(struct relation_term *tu)
{
  memset((char *)tu,0,sizeof(union RelationTermUnion));
}
#endif

static struct {
  long startcheck;
  size_t len;
  size_t cap;
  struct relation_term **buf;
  unsigned long *termstack;
  unsigned long termstackcap;
  long endcheck;
} g_term_ptrs = {1234567890,0,0,NULL,NULL,0,987654321};

#define TPBUF_RESET (g_term_ptrs.len=0)
/* forget about all the terms in the buffer */


/*
 * Now one can ask why a pool and a buffer both? Couldn't one just
 * run a big buffer? Well, yes, but how big? Growing a buffer of
 * complete tokens can cause some system allocators to behave very
 * poorly. Growing a vector of pointers to tokens is much less
 * likely to cause the allocator headaches.
 *
 * The pool has a good growth mechanism and can handle tokens.
 * Tradeoff: it is slower to copy the final token data into a
 * fixed array from pool pointers than from a buffer monolith.
 */
#define TPBUF_INITSIZE 1000
/* initial token buffer capacity */
#define TPBUF_GROW 1000
/* token buffer growth rate */

#define RP_LEN 5
#if (SIZEOF_VOID_P == 8)
#define RP_WID 41
#else
#define RP_WID 63
#endif
/* retune rpwid if the size of tokens changes dramatically */
#define RP_ELT_SIZE (sizeof(union RelationTermUnion))
#define RP_MORE_ELTS 5
/* Number of slots filled if more elements needed.
   So if the pool grows, it grows by RP_MORE_ELTS*RP_WID elements at a time. */
#define RP_MORE_BARS 508
/* This is the number of pool bar slots to add during expansion.
   not all the slots will be filled immediately. */

/* This function is called at compiler startup time and destroy at shutdown.
   One could also recall these every time there is a delete all types. */
void InitRelInstantiator(void) {
  if (g_term_pool != NULL || g_term_ptrs.buf != NULL) {
    ASC_PANIC("ERROR: InitRelInstantiator called twice.\n");
  }
  g_term_pool =
    pool_create_store(RP_LEN, RP_WID, RP_ELT_SIZE, RP_MORE_ELTS, RP_MORE_BARS);
  if (g_term_pool == NULL) {
    Asc_Panic(2, "InitRelInstantiator",
              "ERROR: InitRelInstantiator unable to allocate pool.\n");
  }
  g_term_ptrs.buf = (struct relation_term **)
    ASC_NEW_ARRAY_CLEAR(union RelationTermUnion *,TPBUF_INITSIZE);
  /* don't let the above cast fool you about what's in the array */
  if (g_term_ptrs.buf == NULL) {
    Asc_Panic(2, "InitRelInstantiator",
              "ERROR: InitRelInstantiator unable to allocate memory.\n");
  }
  g_term_ptrs.len = 0;
  g_term_ptrs.cap = TPBUF_INITSIZE;
  g_term_ptrs.termstackcap = 200;
  g_term_ptrs.termstack = ASC_NEW_ARRAY(unsigned long,200);
  if (g_term_ptrs.termstack == NULL) {
    Asc_Panic(2, "InitRelInstantiator",
              "ERROR: InitRelInstantiator unable to allocate memory.\n");
  }
}

/* this function returns NULL when newcap is 0 or when
 * it is unable to allocate the space requested.
 */
static unsigned long *realloc_term_stack(unsigned long newcap){
  if (!newcap) {
    if (g_term_ptrs.termstackcap !=0) {
      ascfree(g_term_ptrs.termstack);
      g_term_ptrs.termstack = NULL;
      g_term_ptrs.termstackcap = 0;
    }
  } else {
    if (newcap >= g_term_ptrs.termstackcap) { /*less than means currently ok */
      unsigned long *newbuf;
      newbuf = (unsigned long *)
        ascrealloc(g_term_ptrs.termstack,(sizeof(unsigned long)*newcap));
      if (newbuf!=NULL) {
        g_term_ptrs.termstack = newbuf;
        g_term_ptrs.termstackcap = newcap;
      } else {
    FPRINTF(ASCERR,"Insufficient memory in relation processor\n");
        return NULL;
      }
    }
  }
  return g_term_ptrs.termstack;
}

void DestroyRelInstantiator(void) {
  assert(g_term_ptrs.buf!=NULL);
  assert(g_term_pool!=NULL);
  ascfree(g_term_ptrs.buf);
  g_term_ptrs.buf = NULL;
  g_term_ptrs.cap = g_term_ptrs.len = (size_t)0;
  if (g_term_ptrs.termstackcap != 0) {
    ascfree(g_term_ptrs.termstack);
    g_term_ptrs.termstack = NULL;
    g_term_ptrs.termstackcap = 0;
  }
  pool_destroy_store(g_term_pool);
  g_term_pool = NULL;
}

void ReportRelInstantiator(FILE *f)
{
  assert(g_term_pool!=NULL);
  FPRINTF(f,"RelInstantiator ");
  pool_print_store(f,g_term_pool,0);
  FPRINTF(f,"RelInstantiator buffer capacity: %lu\n",
    (unsigned long)g_term_ptrs.cap);
}

/* The slower expansion process. */
static void ExpandTermBuf(struct relation_term *t) {
  struct relation_term **newbuf;
  newbuf = (struct relation_term **)ascrealloc(g_term_ptrs.buf,
      (sizeof(struct relation_term *)*(g_term_ptrs.cap+TPBUF_GROW)));
  if (newbuf!=NULL) {
    g_term_ptrs.buf = newbuf;
    g_term_ptrs.cap += TPBUF_GROW;
    g_term_ptrs.buf[g_term_ptrs.len] = t;
    g_term_ptrs.len++;
  } else {
    FPRINTF(ASCERR,
            "ERROR: Relation Instantiator unable to allocate memory.\n");
    /* we have ignored the term pointer, but somebody else still has it: pool*/
  }
  return;
}

/* Appends term to buffer. if buffer full and can't expand, forgets term.*/
static void AppendTermBuf(struct relation_term *t) {
  if (g_term_ptrs.len < g_term_ptrs.cap) {
    g_term_ptrs.buf[g_term_ptrs.len++] = t;
  } else {
    ExpandTermBuf(t);
  }
  return;
}

/*------------------------------------------------------------------------------
  FUNCS TO SIMPLIFY POSTFIX TOKEN LIST

  ...before final creation of the token relation array.
*/

/* returns 1 if term is e_zero, e_real=0.0, or e_int=0 */
static int SimplifyTBIsZero(struct relation_term *arg)
{
  if (RelationTermType(arg)==e_real && R_TERM(arg)->value == 0.0) return 1;
  if (RelationTermType(arg)==e_int && I_TERM(arg)->ivalue == 0) return 1;
  if (RelationTermType(arg)==e_zero) return 1;
  return 0;
}

#ifdef THIS_IS_AN_UNUSED_FUNCTION
/* check a termtype, t, for scalarness. return 1 if so, 0 otherwise. */
static int SimplifyTBIsScalar(enum Expr_enum t)
{
  return (t <= TOK_SCALAR_HIGH && t >= TOK_SCALAR_LOW);
}
#endif  /* THIS_IS_AN_UNUSED_FUNCTION */


/* check a termtype, t, for constantness, return 1 if so, 0 otherwise. */
static int SimplifyTBIsConstant(enum Expr_enum t)
{
  return (t <= TOK_CONSTANT_HIGH && t >= TOK_CONSTANT_LOW);
}

#define ZEROTERM(rtp) SimplifyTBIsZero(rtp)
/* check a term pointer, rtp, for scalarness */
#define SCALARTERM(t) SimplifyTBIsScalar(t)
/* check a termtype, t, for scalarness */
#define CONSTANTTERM(t) SimplifyTBIsConstant(t)
/* check a termtype, t, for constantness */

/*
 * Attempt to simplify unary functions.
 * Returns 1 if arg is not constant.
 * Returns 0 if succeeded, in which case *fn is now morphed to a constant term.
 * Returns -1 if arg value/dimens are inconsistent with function fn.
 * Constant arg with numeric value 0 and wild/no dim are coerced quietly
 * where applicable.
 *
 * Cost: O(1).
 */
static int SimplifyTermBuf_Func(struct relation_term *arg,
                struct relation_term *fn)
{
  CONST dim_type *newdim=NULL;
  double rval;
  /* zero constants */
  if (ZEROTERM(arg)) {
    switch(FuncId(F_TERM(fn)->fptr)) {
    case F_LN:
    case F_LOG10:
    case F_ARCCOSH:
      /* illegal argument. caller will whine. */
      return -1;
    case F_EXP:
    case F_COSH:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless())) {
        arg->t = e_nop;
        fn->t = e_int;
        I_TERM(fn)->ivalue = 1;
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    case F_COS:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),TrigDimension())) {
        arg->t = e_nop;
        fn->t = e_int;
        I_TERM(fn)->ivalue = 1;
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    case F_SIN:
    case F_TAN:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),TrigDimension())) {
        arg->t = e_nop;
        fn->t = e_int;
        I_TERM(fn)->ivalue = 0;
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
#ifdef HAVE_ERF
    case F_ERF:
#endif
    case F_SINH:
    case F_ARCSINH:
    case F_TANH:
    case F_ARCTANH:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless())) {
        arg->t = e_nop;
        fn->t = e_int;
        I_TERM(fn)->ivalue = 0; /* dimensionless integer 0 */
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    case F_CUBE:
    {
      newdim = CubeDimension(TermDimensions(arg),1);
      if (newdim != NULL) {
        arg->t = e_nop;
        fn->t = e_real;
        R_TERM(fn)->value = 0.0;
        R_TERM(fn)->dimensions = newdim;
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    }
    case F_CBRT:
    {
      newdim = ThirdDimension(TermDimensions(arg),1);
      if (newdim != NULL) {
        arg->t = e_nop;
        fn->t = e_real;
        R_TERM(fn)->value = 0.0;
        R_TERM(fn)->dimensions = newdim;
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    }
    case F_SQR:
    {
      newdim = SquareDimension(TermDimensions(arg),1);
      if (newdim != NULL) {
        arg->t = e_nop;
        fn->t = e_real;
        R_TERM(fn)->value = 0.0;
        R_TERM(fn)->dimensions = newdim;
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    }
    case F_SQRT:
    {
      newdim = HalfDimension(TermDimensions(arg),1);
      if (newdim != NULL) {
        arg->t = e_nop;
        fn->t = e_real;
        R_TERM(fn)->value = 0.0;
        R_TERM(fn)->dimensions = newdim;
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    }
    case F_ARCSIN:
    case F_ARCTAN:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless())) {
        arg->t = e_nop;
        fn->t = e_real;
        R_TERM(fn)->value = 0.0;
        R_TERM(fn)->dimensions = TrigDimension();
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    case F_ARCCOS:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless())) {
        arg->t = e_nop;
        fn->t = e_real;
        R_TERM(fn)->value = F_PI_HALF;
        R_TERM(fn)->dimensions = TrigDimension();
        return 0;
      } else {
        return -1; /* dimensional incompatibility */
      }
    case F_ABS:
    case F_HOLD:
    {
      newdim = TermDimensions(arg);
      if (newdim != NULL) {
        arg->t = e_nop;
        fn->t = e_real;
        R_TERM(fn)->value = 0.0;
        R_TERM(fn)->dimensions = newdim;
        return 0;
      } else {
        return -1; /* dimensional insanity */
      }
    }
    case F_LNM:
      return 1; /* user could change lnm epsilon. can't simplify. */
    default:
      FPRINTF(ASCERR,"Unrecognized function in relation.\n");
      return 1;
    }
  }
  /* nonzero int or real */
  if( (arg->t == e_int) || (arg->t == e_real) ) {
    newdim = NULL;
    if (arg->t == e_int) {
      rval = (double)I_TERM(arg)->ivalue;
    } else {
      rval = R_TERM(arg)->value;
    }
    switch(FuncId(F_TERM(fn)->fptr)) {
    /* things that take any trig arg, return dimensionless */
    case F_SIN:
    case F_COS:
    case F_TAN:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),TrigDimension())) {
        newdim = Dimensionless();
      } else {
        return -1; /* dimensional incompatibility */
      }
      break; /* go to fixup */
    /* things that require arg >= 1, return dimless */
    case F_ARCCOSH:
      if( rval < 1.0 ) return -1;
      /* fall through */
    case F_LN:
    case F_LOG10:
      if( rval < 0.0 ) return -1;
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless())) {
        newdim = Dimensionless();
      } else {
        return -1; /* dimensional incompatibility */
      }
      break; /* go to fixup */
    /* things that take any exponentiable  arg, return dimensionless */
    case F_EXP:
    case F_SINH:
    case F_COSH:
      if (fabs(rval) > F_LIM_EXP) return -1;
      /* fall through */
    /* things that take any arg, return dimensionless */
    case F_ARCSINH:
    case F_TANH:
#ifdef HAVE_ERG
    case F_ERF:
#endif
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless())) {
        newdim = Dimensionless();
      } else {
        return -1; /* dimensional incompatibility */
      }
      break;
    case F_ARCTANH:
    /* things that take any arg abs <1, return dimensionless */
      if (fabs(rval) < 1.0 && (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless()))) {
        newdim = Dimensionless();
      } else {
        return -1; /* dimensional incompatibility   or range */
      }
      break;
    case F_CUBE:
    {
      newdim = CubeDimension(TermDimensions(arg),1);
      if (newdim == NULL || fabs(rval) > F_LIM_CUBE) {
        return -1; /* dimensional incompatibility */
      }
    }
      break;
    case F_CBRT:
    {
      newdim = ThirdDimension(TermDimensions(arg),1);
      if (newdim == NULL) {
        return -1; /* dimensional incompatibility , range*/
      }
      break;
    }
    case F_SQR:
    {
      newdim = SquareDimension(TermDimensions(arg),1);
      if (newdim == NULL || fabs(rval) > F_LIM_SQR) {
        return -1; /* dimensional incompatibility , range*/
      }
      break;
    }
    case F_SQRT:
    {
      newdim = HalfDimension(TermDimensions(arg),1);
      if (newdim == NULL || rval < 0.0) {
        return -1; /* dimensional incompatibility or range */
      }
      break;
    }
    /* things that take any trig arg, return dimensionless */
    case F_ARCSIN:
    case F_ARCCOS:
      if ( fabs(rval) <= 1.0 && (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless()))) {
        newdim = TrigDimension();
        break;
      } else {
        return -1; /* dimensional incompatibility */
      }
    case F_ARCTAN:
      if (IsWild(TermDimensions(arg)) ||
          SameDimen(TermDimensions(arg),Dimensionless())) {
        newdim = TrigDimension();
        break;
      } else {
        return -1; /* dimensional incompatibility */
      }
    case F_ABS:
    case F_HOLD:
      newdim = TermDimensions(arg);
      break;
    case F_LNM:
      return 1; /* user could change lnm epsilon. can't simplify. */
    default:
      FPRINTF(ASCERR,"Unrecognized function in relation.\n");
      return 1;
    }
    rval = FuncEval(TermFunc(A_TERM(fn)),rval);
    if (floor(rval)==ceil(rval) && SameDimen(newdim,Dimensionless()) &&
        fabs(rval) < MAXINTREAL) {
      fn->t = e_int;
      I_TERM(fn)->ivalue = (long)floor(rval);
    } else {
      fn->t = e_real;
      R_TERM(fn)->value = rval;
      R_TERM(fn)->dimensions = newdim;
    }
    return 0;
  }
  return 1;
}

static int ArgsForToken(enum Expr_enum t) {
  switch (t) {
  case e_nop:
  case e_undefined:
  case e_glassbox:
  case e_blackbox:
  case e_opcode:
  case e_token:
  case e_zero:
  case e_real:
  case e_int:
  case e_var:
    return 0;
  case e_uminus:
  case e_func:
    return 1;
  case e_plus:
  case e_minus:
  case e_times:
  case e_divide:
  case e_power:
  case e_ipower:
  case e_notequal:
  case e_equal:
  case e_less:
  case e_greater:
  case e_lesseq:
  case e_greatereq:
    return 2;
  case e_maximize:
  case e_minimize:
    return 1;
  default:
    FPRINTF(ASCERR,"ArgsForToken called with illegal token type.\n");
    return -1;
  }
}

/**
 * first = SimplifyTermBuf_SubExprLimit(ts,b,start,tt)
 * unsigned long CONST *ts;      current term stack
 * struct relation_term ** CONST b;  global term ptr array
 * unsigned long start;  starting index IN STACK ts to find needed args
 * enum Expr_enum tt;    term type of operator you want the subexpr for
 * long int first;   term stack position of rightmost arg outside subexpr
 *
 * A little function to find the extent of a postfix subexpression for
 * the args of an operator term in the termstack/termbuf processing.
 * Returns -2 if insanity detected. handles nonoperator tt gracefully (-2).
 *
 * e.g. cos(v1+v2) * v3
 * tt = e_times, ts =>
 *  | V1 | V2 | + | cos | V3 | * |
 *                        ^--------start = 3
 * ^--------first = -1
 *
 * e.g. v1 * (v2 + v3)
 * tt = e_plus, ts =>
 *  | V1 | V2 | V3 | + | * |
 *               ^--------start = 2
 *    ^--------first = 0
 *
 * O(n) n= subexpr length.
 */
static long
SimplifyTermBuf_SubExprLimit(unsigned long CONST *ts,
            struct relation_term ** CONST buf,
            unsigned long start,
            enum Expr_enum tt)
{
  long int first, req_args;

  first = start;
  req_args = ArgsForToken(tt);
  if (first < 0) {
    FPRINTF(ASCERR,"SimplifyTermBuf_SubExpr given malformed subexpression.\n");
  }

  while (first >= 0 && req_args >0) {
    switch(buf[ts[first]]->t) {
    case e_zero:
    case e_real:
    case e_int:
    case e_var:
      req_args--;
      break;
    case e_plus:
    case e_minus:
    case e_times:
    case e_divide:
    case e_power:
    case e_ipower:
      req_args++;
      break;
    case e_func:
    case e_uminus:
      break;
    default:
      FPRINTF(ASCERR,
        "SimplifyTermBuf_SubExpr found illegal argument type (%d).\n",
        buf[ts[first]]->t);
      return -2;
    }
    first--;
  }
  if (first < -1) {
    FPRINTF(ASCERR,"SimplifyTermBuf_SubExpr found malformed subexpression.\n");
  }
  return first;
}

#ifndef NDEBUG
/* some functions to keep assert happy when simplification is in debug */
static int check_gt0(unsigned long i) {
  assert(i);
  return 1;
}
static int check_gt1(unsigned long i) {
  assert(i>1);
  return 1;
}
#endif

/**
 * A function to simplify the term buffer before copying it into a
 * postfix array. Only mandatory dim checking is performed.
 * Cost: O(n) where n = blen.
 *
 * This function is rather large, but simply structured, because speed
 * is important.
 * This is postfix simplification on the cheap. It could be more aggressive,
 * but only at potentially quadratic expense.
 *
 * int level;
 * struct relation_term ** CONST b;
 * CONST unsigned long blen;
 * They are the original term buffer array and its starting length.
 * b stays constant, not the data in it!
 *
 * (the following level definitions are mostly vapor. see relation.h for true.
 * level is how far to go in simplification. it is cumulative.
 * level 0 = do nothing.
 * level 1 = constant folding
 * level 2 = zero reductions. A*0 = 0/A =0. A^0=1;
 * level 3 = converting division by constants into multiplication
 * level 4 = distributing constants over simple mult. (V*C2)*C1 --> V*C3
 *
 * As a side effect, any e_power term that can be resolved to having
 * an integer exponent is converted to an e_ipower.
 *
 * This function is designed to simplifications wrt constants that
 * are easy to do in postfix. If you want something more clever, you
 * need to dress up things in infix, simplify, and put back to postfix.
 * Better you than me, bud.
 *
 * At present level > 1 is ignored; we always do 1-3, never 4.
 *
 * All this goes on in the termbuf array leaving null pointers behind.
 * We will compact the array and adjust the length before leaving this
 * function, so you don't have to care about len changing.
 * The termbuf pointers are from the pool, so we do not free them
 * as terms are eliminated.
 *
 * Internal doc:
 * Because C optimizers are pretty damned good, we aren't going to
 * play pointer games, we will just play subscript of b games.
 * Note that in flight we create null pointers in the already
 * visited buffer, but we always have an argument immediately
 * to the left (b[i-1]) of operator b[i]. If b[i] binary, its
 * right arg is b[i-1] and its left arg is the first nonnull
 * entry b[j] to the left of b[i-1] (j<i-1).
 *
 * The buffer is in postfix. We have no infix to maintain yet.
 * Abbreviations in comments:
 *  U - unary operator
 *  B - binary operator
 *  P - any operator
 *  V - e_var arg
 *  A - any arg
 *  C - any constant arg (e_int, e_real)
 *  R - e_real arg
 *  I - e_int arg
 *  N - null pointer
 * While in flight:
 | A | A | A | A | A | A | A |  termbuf
 *                     ^------- top = rightmost we've considered (current).
 | S | S | S | 0 |
 *             ^----next = next free location to put an index in termstack
 */
static unsigned long SimplifyTermBuf(int level,
                 register struct relation_term ** CONST b,
                 CONST unsigned long blen)
{
  register unsigned long next;
  register unsigned long *ts; /* term stack, should we need it */
  unsigned long top;
  long last;
  unsigned long right;
  int early = 0, err;
  CONST dim_type *newdim;
  long ival;
  double rval;

  if ( level < 1 || !blen ) {
    realloc_term_stack(0);
    return blen;
  }
  ts = realloc_term_stack(blen);
  /* stack gets used a lot, so make him locally managed, reusable mem */
  if (ts==NULL) return blen;
  /* at any trip through this loop we must be able to guarantee
   * some simple change, or that the buffer is suitable for
   * cleanup and return, so that we can handle the rogue operators,
   * args cleanly.
   */
  /* check that stack doesn't start with operator */
  /* should check that stack doesn't start pos 1 with binary operator */
  switch (b[0]->t) {
  case e_var:
  case e_int:
  case e_real:
  case e_zero:
    break;
  default:
    FPRINTF(ASCERR,"Compiler cannot simplify malformed expression\n");
    return blen;
  }

#ifdef NDEBUG
# define TS_TOP (ts[next-1]) /* term address last pushed */
# define TS_LEFT (ts[next-2])
   /* left hand term address IFF current term is binary and the term at TS_TOP is scalar (not operator) */
# define TS_SHIFTPOP ts[next-2] = ts[next-1]; next-- /* overwrite ts_left with ts_top and pop */
#else
# define TS_TOP (check_gt0(next),ts[next-1]) /* term address last pushed */
# define TS_LEFT (check_gt1(next),ts[next-2]) /* left hand term address IFF current term is binary and the term at TS_TOP is scalar (not operator) */
# define TS_SHIFTPOP assert(next>1); ts[next-2] = ts[next-1]; next-- /* overwrite ts_left with ts_top and pop */
#endif
/* keep the above definitions in sync. only difference should be assert. */

#define TS_PUSH(index) ts[next]=(index); next++ /* add a term to the stack */
#define TS_POP next-- /* backup the stack */
#define TS_POP2 next -= 2 /* backup the stack 2 spots */

  for (next=top=0; top < blen; top++) {
    /* pass through the tokens pointers array */
    if (b[top]==NULL) continue; /* so we can go through again if we like */
    /* each case and nested case should be complete in itself for
       readability. do not use fall throughs */
    switch (b[top]->t) {
    case e_var:
    case e_int:
    case e_real:
    case e_zero:
      TS_PUSH(top);
      break;
    case e_nop:
      b[top] = NULL; /* forget nop */
      break;
    case e_func:
      if ( CONSTANTTERM(b[TS_TOP]->t) ) {
    /* C U -> C' */
        if ( (err = SimplifyTermBuf_Func(b[TS_TOP],b[top]) ) != 0 ) {
          /* not simplified. just push later. whine if needed. */
          if (err < 0) {
            FPRINTF(ASCERR,
              "Can't simplify inconsistent argument to unary function.\n");
          }
        } else {
          b[TS_TOP] = NULL; /* kill old arg, func term was morphed. */
          TS_POP;   /* set up to push morphed func in place of arg */
        }
      }
      TS_PUSH(top); /* for all cases in the ifs */
      break;
    case e_uminus:
      switch (b[TS_TOP]->t) {
      case e_int:
        I_TERM(b[TS_TOP])->ivalue = -(I_TERM(b[TS_TOP])->ivalue);
        b[top] = b[TS_TOP]; /* I - => -I */
        b[TS_TOP] = NULL;
        TS_POP;
        TS_PUSH(top);
        break;
      case e_real:
        R_TERM(b[TS_TOP])->value = -(R_TERM(b[TS_TOP])->value);
        b[top] = b[TS_TOP]; /* R - => -R */
        b[TS_TOP] = NULL;
        TS_POP;
        TS_PUSH(top);
        break;
      case e_zero:
        b[top] = b[TS_TOP]; /* -0 = 0 */
        b[TS_TOP] = NULL;
        TS_POP;
        TS_PUSH(top);
        break;
      default:  /* aren't going to distribute or do other funky things */
        TS_PUSH(top);
        break;
      }
      break;

    case e_plus:
      /* A 0 + => NULL NULL A */
      if ( ZEROTERM(b[TS_TOP]) ) {
        /*
         * Note: we really should be checking the dimens of A to match
         * with dimens of 0 if e_real, but we are can't yet.
         */
        b[top] = b[TS_LEFT];    /* overwrite the + with the A */
        b[TS_LEFT] = NULL;  /* null the A old location */
        b[TS_TOP] = NULL;   /* null old location of 0 */
        TS_POP2;
        TS_PUSH(top);
        break;
      }
      switch (b[TS_TOP]->t) {
      case e_var:
        if ( ZEROTERM(b[TS_LEFT]) ) {
          /* 0 V + => NULL NULL V */
          /*
           * Note: we really should be checking the dimens of V to match
           * with dimens of 0 if e_real, but we are don't yet.
           */
          b[TS_LEFT] = NULL; /* null the zero term */
          b[top] = b[TS_TOP]; /* overwrite the + with the V */
          b[TS_TOP] = NULL; /* null old location of V */
          TS_POP2;
          TS_PUSH(top);
        } else {
          TS_PUSH(top);
        }
        break;
      /* 2 constant args? mangle C1 C2 + => C3 of appropriate type,if ok. */
      case e_int:  /* 0 I +, R I +, I I + */
        if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
          /* 2 constant args. mangle C2 I1 + => C3 of appropriate type,if ok.*/
          if (b[TS_LEFT]->t==e_zero) { /* 0 I + */
            b[top] = b[TS_TOP]; /* overwrite the + with the I */
            b[TS_LEFT] = NULL;  /* null the 0 old location */
            b[TS_TOP] = NULL;   /* null old location of I */
            TS_POP2;
            TS_PUSH(top);
            break;
          }
          if (b[TS_LEFT]->t == e_int) { /* I2 I1 + */
            I_TERM(b[TS_TOP])->ivalue += I_TERM(b[TS_LEFT])->ivalue;
            b[top] = b[TS_TOP];   /* overwrite the + with the I1' */
            b[TS_LEFT] = NULL;    /* null the I2 old location */
            b[TS_TOP] = NULL;     /* null old location of I1 */
            TS_POP2;
            TS_PUSH(top);
            break;
          }
          if ( b[TS_LEFT]->t==e_real &&
               ( SameDimen(R_TERM(b[TS_LEFT])->dimensions,Dimensionless()) ||
                 (IsWild(R_TERM(b[TS_LEFT])->dimensions) &&
                  R_TERM(b[TS_LEFT])->value == 0.0)
               )
             ) { /* Ri2(possibly wild 0.0)  I1 + => I1' */
            if (floor(R_TERM(b[TS_LEFT])->value) ==
                  ceil(R_TERM(b[TS_LEFT])->value) &&
                  fabs(R_TERM(b[TS_LEFT])->value) < MAXINTREAL) {
              I_TERM(b[TS_TOP])->ivalue +=
                (long)floor(R_TERM(b[TS_LEFT])->value);
              b[top] = b[TS_TOP];   /* overwrite the + with the I1' */
              b[TS_LEFT] = NULL;    /* null the R2 old location */
              b[TS_TOP] = NULL;     /* null old location of I1 */
              TS_POP2;
              TS_PUSH(top);
              break;
            } else { /* morph + to R' */
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions = Dimensionless();
              R_TERM(b[top])->value =
                R_TERM(b[TS_LEFT])->value + (double)I_TERM(b[TS_TOP])->ivalue;
              b[TS_LEFT] = NULL;    /* null the R2 old location */
              b[TS_TOP] = NULL;     /* null old location of I1 */
              TS_POP2;
              TS_PUSH(top);
              break;
            }
          } else {  /* dimensional conflict can't be fixed */
            FPRINTF(ASCERR,
              "Can't simplify dimensionally inconsistent arguments to +.\n");
            TS_PUSH(top);
          }
          break;
        } else {    /* non C TS_LEFT */
          TS_PUSH(top);
        }
        break;
      case e_real: /* 0 R +, R R +, I R + */
        if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
          /* 2 constant args. mangle C2 R1 + => C3 of appropriate type,if ok.*/
          newdim = CheckDimensionsMatch(TermDimensions(b[TS_TOP]),
                                        TermDimensions(b[TS_LEFT]));
          if (newdim == NULL) {
            FPRINTF(ASCERR,
              "Can't simplify dimensionally inconsistent arguments to +.\n");
            TS_PUSH(top);
            break;
          }
          if (b[TS_LEFT]->t==e_zero) { /* 0 R + */
            b[top] = b[TS_TOP]; /* overwrite the + with the R */
            b[TS_LEFT] = NULL;  /* null the 0 old location */
            b[TS_TOP] = NULL;   /* null old location of R */
            TS_POP2;
            TS_PUSH(top);
            /* if R was wild, it stays wild */
            break;
          }
          if (b[TS_LEFT]->t == e_int) { /* I2 R1 + */
            R_TERM(b[TS_TOP])->value += (double)I_TERM(b[TS_LEFT])->ivalue;
            R_TERM(b[TS_TOP])->dimensions = newdim;
            b[top] = b[TS_TOP];   /* overwrite the + with the R1' */
            b[TS_LEFT] = NULL;    /* null the I2 old location */
            b[TS_TOP] = NULL;     /* null old location of R1 */
            TS_POP2;
            TS_PUSH(top);
            /* if R wild, R becomes dimensionless */
            break;
          }
          if ( b[TS_LEFT]->t==e_real ) { /* R2 R1 + => R1', if R1' whole->I1'*/
            b[top]->t = e_real;
            R_TERM(b[top])->dimensions = newdim;
            R_TERM(b[top])->value =
              R_TERM(b[TS_LEFT])->value + R_TERM(b[TS_TOP])->value;
            b[TS_LEFT] = NULL;    /* null the R2 old location */
            b[TS_TOP] = NULL;     /* null old location of R1 */
            TS_POP2;
            TS_PUSH(top);
            /* if integer valued dimless real, convert to int */
            if (floor(R_TERM(b[top])->value) == ceil(R_TERM(b[top])->value)
                 && SameDimen(R_TERM(b[top])->dimensions,Dimensionless()) &&
                 fabs(R_TERM(b[top])->value) < MAXINTREAL) {
              I_TERM(b[top])->ivalue = (long)R_TERM(b[top])->value;
              b[top]->t = e_int;
            }
            break;
          } else {  /* dimensional conflict can't be fixed */
            FPRINTF(ASCERR,
              "Can't simplify dimensionally inconsistent arguments to +.\n");
            TS_PUSH(top);
          }
          break;
        } else {    /* non C TS_LEFT */
          TS_PUSH(top);
        }
        break; /* end eplus, right arg is e_real */
      default: /* tstop is not 0, R, I, V */
        TS_PUSH(top);
        break;
      } /* end argtype switch of e_plus */
      break;

    case e_minus:
      /* A 0 - => NULL NULL A */
      if ( ZEROTERM(b[TS_TOP]) ) {
        /*
         * Note: we really should be checking the dimens of A to match
         * with dimens of 0 if e_real, but we are can't yet.
         */
        b[top] = b[TS_LEFT];    /* overwrite the - with the A */
        b[TS_LEFT] = NULL;  /* null the A old location */
        b[TS_TOP] = NULL;   /* null old location of 0 */
        TS_POP2;
        TS_PUSH(top);
        break;
      }
      switch (b[TS_TOP]->t) {
      case e_var:
        if ( ZEROTERM(b[TS_LEFT]) ) {
          /* 0 V - => NULL V uminus */
        /*
         * Note: we really should be checking the dimens of V to match
         * with dimens of 0 if e_real, but we are don't yet.
         */
          b[TS_LEFT] = NULL;    /* null the zero term */
          b[top]->t = e_uminus; /* morph - to uminus */
          TS_SHIFTPOP;      /* reduce 0 V => V */
          TS_PUSH(top);
        } else {
          TS_PUSH(top);
        }
        break;
      /* 2 constant args? mangle C1 C2 - => C3 of appropriate type,if ok. */
      case e_int:  /* 0 I -, R I -, I I - */
        if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
          /* 2 constant args. mangle C2 I1 - => C3 of appropriate type,if ok.*/
          if (b[TS_LEFT]->t==e_zero) { /* 0 I - */
            b[top] = b[TS_TOP]; /* overwrite the - with -I */
            I_TERM(b[top])->ivalue = -(I_TERM(b[top])->ivalue);
            b[TS_LEFT] = NULL;  /* null the 0 old location */
            b[TS_TOP] = NULL;   /* null old location of I */
            TS_POP2;
            TS_PUSH(top);
            break;
          }
          if (b[TS_LEFT]->t == e_int) { /* I2 I1 - */
            I_TERM(b[TS_TOP])->ivalue =
              I_TERM(b[TS_LEFT])->ivalue - I_TERM(b[TS_TOP])->ivalue;
            b[top] = b[TS_TOP];   /* overwrite the - with the I1' */
            b[TS_LEFT] = NULL;    /* null the I2 old location */
            b[TS_TOP] = NULL;     /* null old location of I1 */
            TS_POP2;
            TS_PUSH(top);
            break;
          }
          if ( b[TS_LEFT]->t==e_real &&
               ( SameDimen(R_TERM(b[TS_LEFT])->dimensions,Dimensionless()) ||
                 (IsWild(R_TERM(b[TS_LEFT])->dimensions) &&
                  R_TERM(b[TS_LEFT])->value == 0.0)
               )
             ) { /* Ri2(possibly wild 0.0)  I1 - => I1' */
            if (floor(R_TERM(b[TS_LEFT])->value) ==
                  ceil(R_TERM(b[TS_LEFT])->value) &&
                fabs(R_TERM(b[TS_LEFT])->value) < MAXINTREAL) {
              I_TERM(b[TS_TOP])->ivalue =
                (long)floor(R_TERM(b[TS_LEFT])->value)
                 - I_TERM(b[TS_TOP])->ivalue;
              b[top] = b[TS_TOP];   /* overwrite the + with the I1' */
              b[TS_LEFT] = NULL;    /* null the R2 old location */
              b[TS_TOP] = NULL;     /* null old location of I1 */
              TS_POP2;
              TS_PUSH(top);
              break;
            } else { /* morph - to R' */
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions = Dimensionless();
              R_TERM(b[top])->value =
                R_TERM(b[TS_LEFT])->value - (double)I_TERM(b[TS_TOP])->ivalue;
              b[TS_LEFT] = NULL;    /* null the R2 old location */
              b[TS_TOP] = NULL;     /* null old location of I1 */
              TS_POP2;
              TS_PUSH(top);
              break;
            }
          } else {  /* dimensional conflict can't be fixed */
            FPRINTF(ASCERR,
              "Can't simplify dimensionally inconsistent arguments to -.\n");
            TS_PUSH(top);
          }
          break;
        } else {    /* non C TS_LEFT */
          TS_PUSH(top);
        }
        break;

      case e_real: /* 0 R -, R R -, I R - */
        if ( CONSTANTTERM(b[TS_LEFT]->t) ) {
          /* 2 constant args. mangle C2 R1 - => C3 of appropriate type,if ok.*/
          newdim = CheckDimensionsMatch(TermDimensions(b[TS_TOP]),
                                        TermDimensions(b[TS_LEFT]));
          if (newdim == NULL) {
            FPRINTF(ASCERR,
              "Can't simplify dimensionally inconsistent arguments to -.\n");
            TS_PUSH(top);
            break;
          }
          if (b[TS_LEFT]->t==e_zero) { /* 0 R - */
            b[top] = b[TS_TOP]; /* overwrite the - with the R */
            R_TERM(b[top])->value = -(R_TERM(b[top])->value);
            b[TS_LEFT] = NULL;  /* null the 0 old location */
            b[TS_TOP] = NULL;   /* null old location of R */
            TS_POP2;
            TS_PUSH(top);
            /* if R was wild, it stays wild */
            break;
          }
          if (b[TS_LEFT]->t == e_int) { /* I2 R1 - */
            R_TERM(b[TS_TOP])->value =
              (double)I_TERM(b[TS_LEFT])->ivalue - R_TERM(b[TS_TOP])->value;
            R_TERM(b[TS_TOP])->dimensions = newdim;
            b[top] = b[TS_TOP];   /* overwrite the - with the R1' */
            b[TS_LEFT] = NULL;    /* null the I2 old location */
            b[TS_TOP] = NULL;     /* null old location of R1 */
            TS_POP2;
            TS_PUSH(top);
            /* if R wild, R becomes dimensionless */
            break;
          }
          if ( b[TS_LEFT]->t==e_real ) { /* R2 R1 - => R1', if R1' whole->I1'*/
            b[top]->t = e_real;     /* morph - to R */
            R_TERM(b[top])->dimensions = newdim;
            R_TERM(b[top])->value =
              R_TERM(b[TS_LEFT])->value - R_TERM(b[TS_TOP])->value;
            b[TS_LEFT] = NULL;    /* null the R2 old location */
            b[TS_TOP] = NULL;     /* null old location of R1 */
            TS_POP2;
            TS_PUSH(top);
            /* if integer valued dimless real, convert to int */
            if (floor(R_TERM(b[top])->value) == ceil(R_TERM(b[top])->value)
                && SameDimen(R_TERM(b[top])->dimensions,Dimensionless())
                && fabs(R_TERM(b[top])->value) < MAXINTREAL) {
              I_TERM(b[top])->ivalue = (long)R_TERM(b[top])->value;
              b[top]->t = e_int;
            }
            break;
          } else {  /* dimensional conflict can't be fixed */
            FPRINTF(ASCERR,
              "Can't simplify dimensionally inconsistent arguments to -.\n");
            TS_PUSH(top);
          }
          break;
        } else {    /* non C TS_LEFT */
          TS_PUSH(top);
        }
        break; /* end eminus, right arg is e_real */
      default: /* tstop is not 0, R, I, V */
        TS_PUSH(top);
        break;
      } /* end argtype switch of e_minus */
      break;

    case e_times:
      /* needs completing. only C*C done at present. need A*0 reductions */
      if ( !CONSTANTTERM(b[TS_LEFT]->t) && !CONSTANTTERM(b[TS_TOP]->t) ) {
        /* no constants in sight, go on fast. */
        TS_PUSH(top);
        break;
      } else {
        /* some constants in sight, try things. */
        if (b[TS_LEFT]->t == e_zero || b[TS_TOP]->t == e_zero) {
          /* end of A 0 * and 0 A *  => 0 */
          ival = SimplifyTermBuf_SubExprLimit(ts,b,next-1,e_times);
          if ( ival > -2 ) {
            for (last = next-1; last > ival; last--) {
              b[ts[last]] = NULL; /* kill the subexpression tokens */
            }
            next = ival + 1; /* big stack pop */
            b[top]->t = e_zero;
            R_TERM(b[top])->dimensions = WildDimension();
            R_TERM(b[top])->value = 0.0;
            TS_PUSH(top);
            break;
          } else {
            /* we had an error in subexpression limit search */
            TS_PUSH(top);
            break;
          }
        } /* end of A 0 * and 0 A * */
        /* NOTE: here we should be watching for 0.0 e_real and 0 e_int,
         * but as yet we don't have the dimen derivation utility to
         * check these cases and derive a properly dimensioned e_real 0.
         * We are not going to do a dimensionally incorrect shortcut
         * implementation. BAA 3/96
         */
        if ( CONSTANTTERM(b[TS_LEFT]->t) ) { /* C A * =?=> ?*/
          /* LEFT is now ereal or e_int because it passed the 0 and C tests */
          if ( b[TS_TOP]->t == e_real) { /* C R * => C */
            if ( b[TS_LEFT]->t == e_real ) { /* R R * => R */
              newdim = SumDimensions(R_TERM(b[TS_TOP])->dimensions,
                                     R_TERM(b[TS_LEFT])->dimensions,1);
              if ( newdim == NULL || IsWild(newdim) ) { /* bad dim */
                FPRINTF(ASCERR,
               "Mult. by wild or fractional dimension constant not folded.\n");
                TS_PUSH(top);
                break;
              } else { /* dim ok. morph etimes to be result. */
                rval = R_TERM(b[TS_TOP])->value * R_TERM(b[TS_LEFT])->value;
                /* god help us if this overflows... */
                b[top]->t = e_real;
                R_TERM(b[top])->dimensions = newdim;
                R_TERM(b[top])->value = rval;
                b[TS_TOP] = NULL;
                b[TS_LEFT] = NULL;
                TS_POP2;
                TS_PUSH(top);
                break;
              }
            } else { /* I R * => R */
              rval =
                R_TERM(b[TS_TOP])->value * (double)I_TERM(b[TS_LEFT])->ivalue;
              /* god help us if this overflows... */
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions = R_TERM(b[TS_TOP])->dimensions;
              R_TERM(b[top])->value = rval;
              b[TS_TOP] = NULL;
              b[TS_LEFT] = NULL;
              TS_POP2;
              TS_PUSH(top);
              break;
            }
#ifndef NDEBUG
            FPRINTF(ASCERR,"Unexpected error in Simplification (1).\n");
            /* NOT REACHED */
            break;
#endif
          }
          if ( b[TS_TOP]->t == e_int) { /* C I * => C */
            if ( b[TS_LEFT]->t == e_real ) { /* R I * => R */
              rval =
                R_TERM(b[TS_LEFT])->value * (double)I_TERM(b[TS_TOP])->ivalue;
              /* god help us if this overflows... */
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions = R_TERM(b[TS_LEFT])->dimensions;
              R_TERM(b[top])->value = rval;
              b[TS_TOP] = NULL;
              b[TS_LEFT] = NULL;
              TS_POP2;
              TS_PUSH(top);
              break;
            } else { /* I I * => I */
              rval = (double)I_TERM(b[TS_TOP])->ivalue *
                     (double)I_TERM(b[TS_LEFT])->ivalue;
              if (fabs(rval) < (double)(LONG_MAX/2)) {/*result safely integer*/
                b[top]->t = e_int;
                I_TERM(b[top])->ivalue =
                  I_TERM(b[TS_TOP])->ivalue * I_TERM(b[TS_LEFT])->ivalue;
                b[TS_TOP] = NULL;
                b[TS_LEFT] = NULL;
                TS_POP2;
                TS_PUSH(top);
                break;
              } else {
                b[top]->t = e_real;
                R_TERM(b[top])->dimensions = Dimensionless();
                R_TERM(b[top])->value = rval;
                b[TS_TOP] = NULL;
                b[TS_LEFT] = NULL;
                TS_POP2;
                TS_PUSH(top);
                break;
              }
            }
#ifndef NDEBUG
            FPRINTF(ASCERR,"Unexpected error in Simplification (2).\n");
            /* NOT REACHED */
            break;
#endif
          }
        } /* end all C A * alternatives */
        /* if here, no simplifications worked,
         * though constants exist.
         */
        TS_PUSH(top);
        break;
      } /* end of case e_times outermost if */
#ifndef NDEBUG
      FPRINTF(ASCERR,"Unexpected error in Simplification (3).\n");
      /* NOT REACHED */
      break;
#endif

    case e_divide: /* note: A1 A2 / postfix => A1/A2 infix */
      /* needs completing only does C/C at present. needs to do 0/A. */
      if ( !CONSTANTTERM(b[TS_LEFT]->t) && !CONSTANTTERM(b[TS_TOP]->t) ) {
        /* no constants in sight, go on fast. */
        TS_PUSH(top);
        break;
      } else {
        /* some constants in sight, try things. */
        if (b[TS_LEFT]->t == e_zero && b[TS_TOP]->t != e_zero) {
          /*  0 A /  => 0 */
          ival = SimplifyTermBuf_SubExprLimit(ts,b,next-1,e_divide);
          if ( ival > -2 ) {
            for (last = next-1; last > ival; last--) {
              b[ts[last]] = NULL; /* kill the subexpression tokens */
            }
            next = ival + 1; /* big stack pop, could be pop2 */
            b[top]->t = e_zero;
            R_TERM(b[top])->dimensions = WildDimension();
            R_TERM(b[top])->value = 0.0;
            TS_PUSH(top);
            break;
          } else {
            /* we had an error in subexpression limit search */
            TS_PUSH(top);
            break;
          }
        } /* end of 0 A / */
        /* NOTE: here we should be watching for 0.0 e_real and 0 e_int,
         * but as yet we don't
         * check these cases and derive a properly dimensioned e_real 0.
         * We are not going to do a dimensionally incorrect shortcut
         * implementation. BAA 3/96
         */
        if ( ZEROTERM(b[TS_TOP]) ) {
          /* trap A/0 out */
          FPRINTF(ASCERR,"Division by constant 0 not simplified.\n");
          top = blen;
          early = 1; /* set flag that we punted. */
          TS_PUSH(top);
          break;
        } /* end of  A/0 out */
        if ( CONSTANTTERM(b[TS_LEFT]->t) ) { /* C A / =?=> ?*/
          /* LEFT is now R or I because it passed the 0 and C tests */
          if ( b[TS_TOP]->t == e_real) { /* C R / => C */
            if ( b[TS_LEFT]->t == e_real ) { /* R R / => R */
              newdim = DiffDimensions(R_TERM(b[TS_LEFT])->dimensions,
                                     R_TERM(b[TS_TOP])->dimensions,1);
              if ( newdim == NULL || IsWild(newdim) ) { /* bad dim */
                FPRINTF(ASCERR,
               "Div. by wild or fractional dimension constant not folded.\n");
                TS_PUSH(top);
                break;
              } else { /* dim ok. morph edivide to be result. */
                rval = R_TERM(b[TS_LEFT])->value / R_TERM(b[TS_TOP])->value;
                /* god help us if this overflows/underflows... */
                b[top]->t = e_real;
                R_TERM(b[top])->dimensions = newdim;
                R_TERM(b[top])->value = rval;
                b[TS_TOP] = NULL;
                b[TS_LEFT] = NULL;
                TS_POP2;
                TS_PUSH(top);
                break;
              }
            } else { /* I R / => R */
              rval =
                ((double)I_TERM(b[TS_LEFT])->ivalue) /R_TERM(b[TS_TOP])->value;
              /* god help us if this overflows... */
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions =
                DiffDimensions(Dimensionless(),
                               R_TERM(b[TS_TOP])->dimensions,0);
                /* diff dimens not checked here because top is dimensionless */
              R_TERM(b[top])->value = rval;
              b[TS_TOP] = NULL;
              b[TS_LEFT] = NULL;
              TS_POP2;
              TS_PUSH(top);
              break;
            }
#ifndef NDEBUG
            FPRINTF(ASCERR,"Unexpected error in Simplification (4).\n");
            /* NOT REACHED */
            break;
#endif
          }
          if ( b[TS_TOP]->t == e_int) { /* C I / => C */
            if ( b[TS_LEFT]->t == e_real ) { /* R I / => R */
              rval =
                R_TERM(b[TS_LEFT])->value / (double)I_TERM(b[TS_TOP])->ivalue;
              /* god help us if this overflows... */
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions = R_TERM(b[TS_LEFT])->dimensions;
              R_TERM(b[top])->value = rval;
              b[TS_TOP] = NULL;
              b[TS_LEFT] = NULL;
              TS_POP2;
              TS_PUSH(top);
              break;
            } else { /* I I / => R! Integer division is NOT allowed */
              rval = (double)I_TERM(b[TS_LEFT])->ivalue;
              rval /= (double)I_TERM(b[TS_TOP])->ivalue;
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions = Dimensionless();
              R_TERM(b[top])->value = rval;
              b[TS_TOP] = NULL;
              b[TS_LEFT] = NULL;
              TS_POP2;
              TS_PUSH(top);
              break;
            }
#ifndef NDEBUG
            FPRINTF(ASCERR,"Unexpected error in Simplification (5).\n");
            /* NOT REACHED */
            break;
#endif
          }
        } /* end all C A / alternatives */
        if ( CONSTANTTERM(b[TS_TOP]->t) ) { /* A C / => A (1/C) * */
          /* we screened out 0 above, so its int or real */
          if (b[TS_TOP]->t == e_real) { /* A R / => A R * */
            rval = 1/R_TERM(b[TS_TOP])->value;
            /* god help us if this overflows... */
            b[top]->t = e_times; /* morph / to * */
            /* flip dimens */
            R_TERM(b[TS_TOP])->dimensions =
              DiffDimensions(Dimensionless(),R_TERM(b[TS_TOP])->dimensions,0);
            /* diff dimens not checked here because top is dimensionless */
            R_TERM(b[TS_TOP])->value = rval; /* flip value */
            TS_PUSH(top);
            break;
          }
          if (b[TS_TOP]->t == e_int) { /* A I / => A I * */
            rval = 1.0/(double)I_TERM(b[TS_TOP])->ivalue;
            /* god help us if this overflows... */
            b[top]->t = e_times; /* morph / to * */
            /* flip dimens */
            b[TS_TOP]->t = e_real; /* morph int to real */
            R_TERM(b[TS_TOP])->value = rval; /* flip value */
            R_TERM(b[TS_TOP])->dimensions = Dimensionless();
            TS_PUSH(top);
            break;
          }
        } /* end of morphing A/C => A*c' */
        /* if here, no simplifications worked,
         * though constants exist.
         */
        TS_PUSH(top);
        break;
      } /* end of case e_divide outermost if */
      /* NOT REACHED */
#ifndef NDEBUG
      FPRINTF(ASCERR,"Unexpected error in Simplification (6).\n");
      break;
#endif
    case e_power: /* DANGER! WILL ROBINSON, DANGER! possible fall through */
      /* exponents must be dimensionless to make any sense */
      if (b[TS_TOP]->t == e_zero || b[TS_TOP]->t == e_int ||
           (b[TS_TOP]->t == e_real &&
             ( SameDimen(R_TERM(b[TS_TOP])->dimensions,Dimensionless()) ||
               IsWild(R_TERM(b[TS_TOP])->dimensions) ) &&
            floor(R_TERM(b[TS_TOP])->value)==ceil(R_TERM(b[TS_TOP])->value) &&
            fabs(R_TERM(b[TS_TOP])->value) < MAXINTREAL)
        ) { /* big if ipowerable */
        if (b[TS_TOP]->t == e_real) { /* morph real to int */
          b[TS_TOP]->t = e_int;
          I_TERM(b[TS_TOP])->ivalue = (long)R_TERM(b[TS_TOP])->value;
        }
        /* e_zero and e_int are appropriate to ipower and need no morph */
        b[top]->t = e_ipower; /* morph to ipower and fall through */
        /* FALL THROUGH! FALL THROUGH! FALL THROUGH! FALL THROUGH! */
        /* we aren't supposed to allow fall, but this is really the
           most perfect place to do power=>ipower conversion.
           Note that very large exponent values may be impossible later. */
      } else {
        /* still need to code C^R case. A^0 promoted to ipow, not here */
        if ( CONSTANTTERM(b[TS_LEFT]->t) && CONSTANTTERM(b[TS_TOP]->t) ) {
          /* C is either 0, int, or real. R is nonintegral (or damn big) real.
             Because R is not integer, C must be positive and dimensionless,
             and also small enough not to overflow: C > 1 =>
             check for pow(DBL_MAX,1/R) > R */
          if ( !SameDimen(R_TERM(b[TS_TOP])->dimensions,Dimensionless()) &&
               !IsWild(R_TERM(b[TS_TOP])->dimensions) ) {
            FPRINTF(ASCERR,"Illegal dimensioned exponent found in power.\n");
            top=blen;
            early = 1; /* set flag that we punted. */
            break;
          }
          if (b[TS_LEFT]->t == e_zero) { /* 0^R, R!=0 => 1 */
            b[top]->t = e_int;
            I_TERM(b[top])->ivalue = 1;
            b[TS_TOP] = NULL;
            b[TS_LEFT] = NULL;
            TS_POP2;
            TS_PUSH(top);
            break;
            /* end of 0^R */
          } else {
            if (b[TS_LEFT]->t == e_real) { /* R^R */
              if ( !SameDimen(R_TERM(b[TS_LEFT])->dimensions,Dimensionless())
                   && !IsWild(R_TERM(b[TS_LEFT])->dimensions) ) {
                /* can happen on very large exponents */
                FPRINTF(ASCERR,
                  "Illegal dimensioned base raised to nonintegral power.\n");
                top = blen;
                early = 1; /* set flag that we punted. */
                break;
              } else { /* R(dimless)^R , err if R ln(C) > ln(DBL_MAX) */
                if (R_TERM(b[TS_LEFT])->value < 0) {
                  /* can happen on very large exponents */
                  FPRINTF(ASCERR,
                    "Illegal negative base raised to nonintegral power.\n");
                  top = blen;
                  early = 1; /* set flag that we punted. */
                  break;
                }
                if (R_TERM(b[TS_LEFT])->value == 0.0) {
                  /* R!=0, 0^R = 1 */
                  b[top]->t = e_int;
                  I_TERM(b[top])->ivalue = 0;
                  b[TS_TOP] = NULL;
                  b[TS_LEFT] = NULL;
                  TS_POP2;
                  TS_PUSH(top);
                  break;
                }
                if ( R_TERM(b[TS_TOP])->value*log(R_TERM(b[TS_LEFT])->value) >
                  F_LIM_EXP) { /* overflow in result. let solver do it */
                  TS_PUSH(top);
                  break;
                } else {
                  b[top]->t = e_real;
                  R_TERM(b[top])->dimensions = Dimensionless();
                  R_TERM(b[top])->value =
                    pow(R_TERM(b[TS_LEFT])->value,R_TERM(b[TS_TOP])->value);
                  b[TS_TOP] = NULL;
                  b[TS_LEFT] = NULL;
                  TS_POP2;
                  TS_PUSH(top);
                  break;
                }
              }
              /* end of R^R */
            } else { /* I^R */
              if (I_TERM(b[TS_LEFT])->ivalue < 0) {
                /* can happen on very large exponents */
                FPRINTF(ASCERR,
                  "Illegal negative base raised to nonintegral power.\n");
                top = blen;
                early = 1; /* set flag that we punted. */
                break;
              }
              if (I_TERM(b[TS_LEFT])->ivalue == 0) {
                /* R!=0, 0^R = 1 */
                b[top]->t = e_int;
                I_TERM(b[top])->ivalue = 0;
                b[TS_TOP] = NULL;
                b[TS_LEFT] = NULL;
                TS_POP2;
                TS_PUSH(top);
                break;
              }
              if ( R_TERM(b[TS_TOP])->value *
                   log((double)I_TERM(b[TS_LEFT])->ivalue) > F_LIM_EXP) {
                /* overflow in result. let solver do it */
                TS_PUSH(top);
                break;
              } else {
                b[top]->t = e_real;
                R_TERM(b[top])->dimensions = Dimensionless();
                R_TERM(b[top])->value =
                  pow((double)I_TERM(b[TS_LEFT])->ivalue,
                      R_TERM(b[TS_TOP])->value);
                b[TS_TOP] = NULL;
                b[TS_LEFT] = NULL;
                TS_POP2;
                TS_PUSH(top);
                break;
              }
              /* end of  I^R  */
            }
            /* end of I,R ^R */
          }
          /* end of 0,I,R ^R */
        } else {
          TS_PUSH(top);
          /* remaining A^A2 where A2 => R or V or expr */
        }
        /* end of not morphing to ipower */
        break;
      }
      /* FALL THROUGH if morphing to ipower test succeeded */

    case e_ipower:
      if ( ZEROTERM(b[TS_TOP]) ) {
        /* A^0 */
        if ( ZEROTERM(b[TS_LEFT]) ) {
          /* 0^0 */
          FPRINTF(ASCERR,"0 raised to 0 power is undefined.\n");
          top=blen;
          early = 1; /* set flag that we punted. */
          break;
        } else {
          /* A^0 => 1 */
          ival = SimplifyTermBuf_SubExprLimit(ts,b,next-1,e_ipower);
          if ( ival > -2 ) {
            for (last = next-1; last > ival; last--) {
              b[ts[last]] = NULL; /* kill the subexpression tokens */
            }
            next = ival + 1; /* big stack pop */
            b[top]->t = e_int;
            I_TERM(b[top])->ivalue = 1;
            TS_PUSH(top);
            break;
          } else {
            /* we had an error */
            TS_PUSH(top);
            break;
          }
        }
      } else { /* A^I, I!=0, A!=0  => R or I as needed */
        if (CONSTANTTERM(b[TS_LEFT]->t)) { /* C^I */
          if (b[TS_LEFT]->t == e_real) { /* R^I */
            /* err if I*ln(R) > ln(DBL_MAX) */
            if ( I_TERM(b[TS_TOP])->ivalue*log(fabs(R_TERM(b[TS_LEFT])->value))
                 > F_LIM_EXP) { /* overflow in result. let solver do it */
              TS_PUSH(top);
              break;
            } else {
              ival = I_TERM(b[TS_TOP])->ivalue;
              newdim = PowDimension(ival,R_TERM(b[TS_LEFT])->dimensions,1);
              if (newdim==NULL) {
                FPRINTF(ASCERR,
                  "Dimensional overflow in exponentiation of constant.\n");
                TS_PUSH(top);
                break;
              }
              b[top]->t = e_real;
              R_TERM(b[top])->dimensions = newdim;
              R_TERM(b[top])->value =
                asc_ipow(R_TERM(b[TS_LEFT])->value,(int)ival);
                /* cast of ival is accurate if newdim was ok */
              b[TS_TOP] = NULL;
              b[TS_LEFT] = NULL;
              TS_POP2;
              TS_PUSH(top);
              break;
            }
            /* end of R^I */
          } else { /* I^I */
            ival = I_TERM(b[TS_TOP])->ivalue;
            if ( ival * log((double)abs(I_TERM(b[TS_LEFT])->ivalue))
                 > F_LIM_EXP) {
              /* overflow in result. let solver do it */
              TS_PUSH(top);
              break;
            }
            if (abs(ival) < INT_MAX) { /* this could be a little better */
              rval = asc_ipow((double)I_TERM(b[TS_LEFT])->ivalue,
                              (int)I_TERM(b[TS_TOP])->ivalue);
              if (fabs(rval) > MAXINTREAL || floor(rval)!=ceil(rval) ) {
                b[top]->t = e_real;
                R_TERM(b[top])->dimensions = Dimensionless();
                R_TERM(b[top])->value = rval;
              } else { /* can be an int safely */
                b[top]->t = e_int;
                I_TERM(b[top])->ivalue = (long)rval;
              }
              b[TS_TOP] = NULL;
              b[TS_LEFT] = NULL;
              TS_POP2;
              TS_PUSH(top);
              break;
            } else {
              /* exponent to damn big */
              TS_PUSH(top);
              break;
            }
            /* end of  I^I  */
          } /* end of C^I */
        } else {
          TS_PUSH(top);
          break;
        }
#ifndef NDEBUG
        FPRINTF(ASCERR,"Unexpected error in Simplification (7).\n");
        break; /* NOT REACHED */
#endif
      }
#ifndef NDEBUG
      FPRINTF(ASCERR,"Unexpected error in Simplification (8).\n");
      break; /* NOT REACHED */
#endif
    /* end e_ipower */

   /* all the following are bogus in instantiated tokens at this time. (2/96)
    * e_subexpr,e_const,e_par,
    * e_card,e_choice,e_sum,e_prod,e_union,e_inter,e_in,e_st,
    * e_glassbox,e_blackbox,e_opcode,e_token,
    * e_or,e_and,e_boolean,e_set,e_symbol,
    * e_equal,e_notequal,e_less,e_greater,e_lesseq,e_greatereq,e_not,
    * e_qstring,
    * e_maximize,e_minimize,
    * e_undefined
    */
    default:
      FPRINTF(ASCERR,"Unexpected token in relation simplification.\n");
      FPRINTF(ASCERR,"Returning early.\n");
      top=blen;
      early = 1; /* flag that we punted. */
      break;
    }
  }
  /*
   *  cleanup reduced expression, if needed.
   *  after the for loop, next is now the length of the postfix expression,
   *  or garbage if early is true.
   */
  if (blen <= next) return blen; /* no simplification, oh well. */
  right = 0;
  while (right < blen && b[right] != NULL) right++; /* find first null */
  for(last = right; right < blen; right++) { /* shift nonnulls left */
    if (b[right] != NULL) {
      b[last] = b[right];
      last++;
    }
  }
  if (!early && last != (long)next) {
    FPRINTF(ASCERR,"Confusing token counts in Simplify\n");
  }
  right = last;
  while (last<(long)blen) { /* null remainder, if any, of pointers */
    b[last] = NULL;
    last++;
  }
  return right;
}
/* END SimplifyTermBuf */


struct relation_side_temp {
  unsigned long length;
  union RelationTermUnion *side;
};

static struct relation_term
*InfixArr_MakeSide(CONST struct relation_side_temp *, int *);
/* forward declaration */

/** returns 1 if converting buf is successful
 * returns 0 if buf empty or insufficient memory.
 * The structure tmp given is filled with an array of terms
 * and its length. You must free the array if you decide you
 * don't want it. We don't care how the structure is initialized.
 */
static int ConvertTermBuf(struct relation_side_temp *tmp)
{
  union RelationTermUnion *arr = NULL;
  unsigned long len,c;

  len = SimplifyTermBuf(g_simplify_relations,g_term_ptrs.buf,g_term_ptrs.len);
  if (len < 1) return 0;
  arr = ASC_NEW_ARRAY(union RelationTermUnion,len);
  if (arr==NULL) {
    FPRINTF(ASCERR,"Create Token Relation: Insufficient memory :-(.\n");
    return 0;
  }
  for (c=0; c<len; c++) {
    arr[c] = *(UNION_TERM(g_term_ptrs.buf[c]));
  }
  tmp->side = arr;
  tmp->length = len;
  return 1;
}

/**
 *  usually we want to reset both simultaneously. reset our
 *  pooling and buffering data.
 */
static
void DestroyTermList(void) {
  POOL_RESET;
  TPBUF_RESET;
}

/**
    create a term from the pool
*/
static struct relation_term *CreateOpTerm(enum Expr_enum t)
{
  struct relation_term *term;
  term = POOL_ALLOCTERM;
  assert(term!=NULL);
  PTINIT(term);
  term->t = t;
  if (t==e_uminus) {
    U_TERM(term)->left = NULL;
  } else {
    B_TERM(term)->left = NULL;
    B_TERM(term)->right = NULL;
  }
  return term;
}

/** create a term from the pool, inserting it
 * in pointer sorted order on g_relation_var_list.
 * Note that this and ModifyTokenRelationPointers are the
 * only places where the sort
 * order of the var list matters.
 * In fact, in most cases we could equally afford
 * linear search and that would give us repeatability
 * across platforms and runs since the vars will be
 * then encountered in a constant order determined
 * by how the user wrote the equation.
 * Needs consideration, especially in light of
 * potential to improve relation sharing.
 * In particular, we could then easily share
 * in a fine-grained manner those relations with
 * only a single index involved and no internal sums/products,
 * such as f[i] = x[i]*Ftot; in[i].f = out[i].f;
 * x = hold(x);
 * which could be pretty darn common forms.
 */
static struct relation_term *CreateVarTerm(CONST struct Instance *i)
{
  struct relation_term *term;
  unsigned long pos;
  if (0 != (pos = gl_search(g_relation_var_list,i,(CmpFunc)CmpP))) {
    /* find var if already on relations var list */
    term = POOL_ALLOCTERM;
    assert(term!=NULL);
    PTINIT(term);
    term->t = e_var;
    V_TERM(term) -> varnum = pos;
  } else {
    /* or add it to the var list */
    gl_append_ptr(g_relation_var_list,(VOIDPTR)i);
    term = POOL_ALLOCTERM;
    assert(term!=NULL);
    PTINIT(term);
    term->t = e_var;
    V_TERM(term) -> varnum = gl_length(g_relation_var_list);
  }
  return term;
}

/** create a term from the pool */
static struct relation_term *CreateIntegerTerm(long int v)
{
  struct relation_term *term;
  term = POOL_ALLOCTERM;
  assert(term!=NULL);
  PTINIT(term);
  term->t = e_int;
  I_TERM(term) -> ivalue = v;
  return term;
}

/** create a term from the pool */
static struct relation_term *CreateRealTerm(double v, CONST dim_type *dim)
{
  struct relation_term *term;
  term = POOL_ALLOCTERM;
  assert(term!=NULL);
  PTINIT(term);
  term->t = e_real;
  R_TERM(term) -> value = v;
  R_TERM(term) -> dimensions = dim;
  return term;
}

/** create a term from the pool. Zero terms look like real, wild zeros */
static struct relation_term *CreateZeroTerm(void)
{
  struct relation_term *term;
  term = POOL_ALLOCTERM;
  assert(term!=NULL);
  PTINIT(term);
  term->t = e_zero;
  R_TERM(term)->value = 0.0;
  R_TERM(term)->dimensions = WildDimension();
  return term;
}

/** create a term from the pool */
static struct relation_term *CreateFuncTerm(CONST struct Func *f)
{
  struct relation_term *term;
  term = POOL_ALLOCTERM;
  assert(term!=NULL);
  PTINIT(term);
  term->t = e_func;
  F_TERM(term) -> fptr = f;
  F_TERM(term) -> left = NULL;
  return term;
}

/** create a term from the pool */
#ifdef  THIS_IS_AN_UNUSED_FUNCTION
static struct relation_term *CreateNaryTerm(CONST struct Func *f)
{
  struct relation_term *term;
  term = POOL_ALLOCTERM;
  assert(term!=NULL);
  PTINIT(term);
  term->t = e_func;
  N_TERM(term)->fptr = f;
  N_TERM(term)->args = NULL;
  return term;
}
#endif  /* THIS_IS_AN_UNUSED_FUNCTION */


/**
    This function creates and *must* create the memory
    for the structure and for the union that the structure
    points to.

    Too much code depends on the pre-existence of a properly initialized union.

    If copyunion is crs_NOUNION, the share ptr is init to NULL and user
    must set refcount,relop after the allocate a UNION or whatever.
    If copyunion is crs_NEWUNION, share ptr is allocated and configured.
*/
struct relation *CreateRelationStructure(enum Expr_enum relop,int copyunion)
{
  struct relation *newrelation;

  newrelation = ASC_NEW(struct relation);
  assert(newrelation!=NULL);
  /* CONSOLE_DEBUG("Created 'struct relation' at %p",newrelation); */

  newrelation->residual = DBL_MAX;
  newrelation->multiplier = DBL_MAX;
  newrelation->nominal = 1.0;
  newrelation->iscond = 0;
  newrelation->vars = NULL;
  newrelation->d =(dim_type *)WildDimension();
  newrelation->externalData = NULL;

  if (copyunion) {
    newrelation->share = ASC_NEW(union RelationUnion);
    assert(newrelation->share!=NULL);
    RelationRefCount(newrelation) = 0;
    RelRelop(newrelation) = relop;
#if TOKENDOMINANT
    RTOKEN(newrelation).lhs_term = NULL;
    RTOKEN(newrelation).rhs_term = NULL;
    RTOKEN(newrelation).lhs = NULL;
    RTOKEN(newrelation).rhs = NULL;
    RTOKEN(newrelation).lhs_len = 0;
    RTOKEN(newrelation).rhs_len = 0;
    RTOKEN(newrelation).btable = 0;
    RTOKEN(newrelation).bindex = 0;
#else
    memset((char *)(newrelation->share),0,sizeof(union RelationUnion));
#endif
  } else {
    newrelation->share = NULL;
  }
  return newrelation;
}


/*------------------------------------------------------------------------------
  EXTERNAL CALL PROCESSING
*/

/** @file "relation.h"
    @note
    A special note on external relations

    External relations behave like relations but they also behave like
    procedures. As such when they are constructed and invoked they expect
    a particular ordering of their variables.

    However there are some operations that can mess up (reduce) the number
    of incident variables on the incident varlist -- ATSing 2 variables in the
    *same* relation will do this. BUT we still need to maintain the number
    of variables in the call to the evaluation routine.

    Consider the following example:
    An glassbox relation is constructed as: test1(x[46,2,8,9] ; 2);
    It *requires* 4 arguements, but its incident var count could be anything
    from 1 <= n <= 4, depending on how many ATS are done.
    The ATS/alias will be done even before we have constructed the relation,
    so we just issue warnings.
*/
struct relation *CreateBlackBoxRelation(struct Instance *relinst
        , struct Instance *subject
        , struct gl_list_t *inputs
        , struct BlackBoxCache * common
        , unsigned long lhsIndex
        , CONST char *context
){
  struct relation *result;
  struct gl_list_t *varlist;
  struct BlackBoxData *bbd;
  struct Instance *var = NULL;
  unsigned long *inputArgs;
  int32 inputsLen;
  int argloc;
  unsigned long c,len,pos, lhsVarNumber;
  unsigned long n_inputs;
  CONST dim_type *d;

  /* CONSOLE_DEBUG("CREATING BLACK BOX RELATION"); */

  n_inputs = gl_length(inputs);
  len = n_inputs + 1; /* an extra for the output variable. */

  /*
    Create the BlackBox relation structure.
    output var always first in varlist.
  */
  bbd = CreateBlackBoxData(common);
  inputsLen = BlackBoxCacheInputsLen(common);
  inputArgs = (unsigned long *)ascmalloc(sizeof(unsigned long) * inputsLen);
  lhsVarNumber = 1;
  varlist = gl_create(len);

  /* add the subject */
  gl_append_ptr(varlist,(VOIDPTR)subject); /* add the subject */
  AddRelation(subject,relinst);

  /* now loop, warning of merges and collecting varlist position
     of each arg into argloc.
  */
  argloc = 0;
  for (c=1; c<=n_inputs; c++) {
    var = (struct Instance *)gl_fetch(inputs,c);
    pos = gl_search(varlist,var,(CmpFunc)CmpP);
    switch (pos) {
    case 0:
      gl_append_ptr(varlist,(VOIDPTR)var);
      AddRelation(var,relinst);
      break;
    case 1:
      ERROR_REPORTER_HERE(ASC_USER_WARNING,"In external relation %s[%d],"
        " output %d and input %d are merged."
        " This will probably destroy your chances of achieving convergence,"
        " unless you are very careful."
        , context, (int)lhsIndex+1, (int)lhsIndex, (int)pos
      );
      break;
    default:
      ERROR_REPORTER_HERE(ASC_USER_WARNING,"In external relation %s[%d],"
        " input %d and input %d are merged."
        " This will probably destroy your chances of achieving convergence,"
        " unless you are very careful."
        , context, (int)lhsIndex+1, argloc, (int)pos
      );
      break;
    }
    if (pos) {
      inputArgs[argloc] = pos;
    } else {
      inputArgs[argloc] = gl_length(varlist);
    }
    argloc++;
  }

  /*
    Now make the main relation structure and put it all
    together.
  */
  result = CreateRelationStructure(e_equal,crs_NEWUNION);
  RelationRefCount(result) = 1;
  result->externalData = bbd;
  result->vars = varlist;
  RBBOX(result).inputArgs = inputArgs;
  RBBOX(result).lhsindex = lhsIndex;
  RBBOX(result).lhsvar = lhsVarNumber;
  d = RealAtomDims(subject);
  if (!IsWild(d)) {
    SetRelationDim(result,d);
  }
  return result;
}


struct relation *CreateGlassBoxRelation(struct Instance *relinst,
                    struct ExternalFunc *efunc,
                    struct gl_list_t *varlist,
                    int gbindex,
                    enum Expr_enum relop)
{
  struct relation *result;
  struct Instance *var;
  struct gl_list_t *newlist = NULL;
  int *tmp = NULL, *args = NULL;
  unsigned long len,c,pos;

  len  = gl_length(varlist);
  /*
   * Make the variables aware that they are incident
   * in this relation instance. At the same time set up
   * the args list indexing.
   */
  if (len) {
    tmp = args = ASC_NEW_ARRAY_CLEAR(int,len+1);
    newlist = gl_create(len);

    for (c=1;c<=len;c++) {
      var = (struct Instance *)gl_fetch(varlist,c);
      pos = gl_search(newlist,var,(CmpFunc)CmpP);
      if (pos) {
    FPRINTF(ASCERR,"Incidence for external relation will be inaccurate\n");
    *tmp++ = (int)pos;
      }
      else{
    gl_append_ptr(newlist,(VOIDPTR)var);
    *tmp++ = (int)gl_length(newlist);
    AddRelation(var,relinst);
      }
    }
  }
  *tmp = 0; /* terminate */

  /*
    Create the relation data structure and append the
    varlist.
  */
  result = CreateRelationStructure(relop,crs_NEWUNION);
  RelationRefCount(result) = 1;
  RGBOX(result).efunc = efunc;
  RGBOX(result).args = args;
  RGBOX(result).nargs = (int)len;
  RGBOX(result).index = gbindex;
  result->vars = newlist;
  return result;
}

/*------------------------------------------------------------------------------
  TOKENRELATION PROCESSING AND GENERAL EXPR-TO-RELATION CHECK ROUTINES
*/

static
struct value_t CheckIntegerCoercion(struct value_t v)
{
  if ((ValueKind(v)==real_value) && (RealValue(v)==0.0) &&
      IsWild(RealValueDimensions(v)) ){
    DestroyValue(&v);
    return CreateIntegerValue(0,1); /* assume this is a constant then */
  }
  else return v;
}

static
int ProcessListRange(CONST struct Instance *ref,
             CONST struct Expr *low,
             CONST struct Expr *up,
             int *added,
             int i,
             enum relation_errors *err,
             enum find_errors *ferr)
{
  struct value_t lower,upper;
  struct relation_term *term;
  long lv,uv;
  assert(GetEvaluationContext()==NULL);
  SetEvaluationContext(ref);
  lower = EvaluateExpr(low,NULL,InstanceEvaluateName);
  upper = EvaluateExpr(up,NULL,InstanceEvaluateName);
  SetEvaluationContext(NULL);
  lower = CheckIntegerCoercion(lower);
  upper = CheckIntegerCoercion(upper);
  if ((ValueKind(lower)==integer_value)&&(ValueKind(upper)==integer_value)){
    lv = IntegerValue(lower);
    uv = IntegerValue(upper);
    while(lv<=uv){
      term = CreateIntegerTerm(lv);
      AppendTermBuf(term);
      if ((*added)++) {
    switch(i){
    case SUM:
      term = CreateOpTerm(e_plus);
      break;
    case PROD:
      term = CreateOpTerm(e_times);
      break;
    }
    AppendTermBuf(term);
      }
      lv++;
    }
    return 0;
  }
  else{
    if(ValueKind(lower)==error_value) {
      FigureOutError(lower,err,ferr);
      return 1;
    }
    if(ValueKind(upper)==error_value){
      FigureOutError(upper,err,ferr);
      return 1;
    }
    *err = incorrect_structure;
    FPRINTF(ASCERR,"incorrect_structure in ProcessListRange\n");
    return 1;
  }
}

static
CONST struct Expr *ExprContainsSuchThat(register CONST struct Expr *ex)
{
  while(ex!=NULL){
    if (ExprType(ex)==e_st) return ex;
    ex = NextExpr(ex);
  }
  return ex;
}

/**
 *  Here we give up if vars are not well defined.
 *  At present e_var acceptable ARE:
 *  REAL_ATOM_INSTANCE
 *  Well defined Real and Integer constants.
 *  Everything else is trash.
 *  CreateTermFromInst() and CheckExpr() must have matching semantics.
 */
static
struct relation_term *CreateTermFromInst(struct Instance *inst,
                     struct Instance *rel,
                     enum relation_errors *err)
{
  struct relation_term *term;
  switch(InstanceKind(inst)){
  case REAL_ATOM_INST:
    term = CreateVarTerm(inst);
    AddRelation(inst,rel);
    return term;
  case REAL_CONSTANT_INST:
    if ( AtomAssigned(inst) && !IsWild(RealAtomDims(inst)) ){
      term = CreateRealTerm(RealAtomValue(inst),RealAtomDims(inst));
      return term;
    }
    else{
      if ( IsWild(RealAtomDims(inst)) && AtomAssigned(inst) ) {
    *err = real_value_wild;
      } else {
    *err = real_value_undefined;
      }
      return NULL;
    }
  case INTEGER_CONSTANT_INST:
    if (AtomAssigned(inst)){
      term = CreateIntegerTerm(GetIntegerAtomValue(inst));
      return term;
    }
    else{
      *err = integer_value_undefined;
      return NULL;
    }
  case REAL_INST:
    *err = incorrect_real_inst_type;
    return NULL;
  case INTEGER_ATOM_INST:
  case INTEGER_INST:
    *err = incorrect_integer_inst_type;
    return NULL;
  case SYMBOL_ATOM_INST:
  case SYMBOL_CONSTANT_INST:
  case SYMBOL_INST:
    *err = incorrect_symbol_inst_type;
    return NULL;
  case BOOLEAN_ATOM_INST:
  case BOOLEAN_CONSTANT_INST:
  case BOOLEAN_INST:
    *err = incorrect_boolean_inst_type;
    return NULL;
  default:
    *err = incorrect_inst_type;
    return NULL;
  }
}

/* forward declaration */
static int AppendList( CONST struct Instance *,
         struct Instance *,
         CONST struct Set *,
         int ,
         enum relation_errors *,
         enum find_errors *);

/**
    @todo document this

    Convert a part of an expression into part of a relation (in postfix)?
*/
static
int ConvertSubExpr(CONST struct Expr *ptr,
           CONST struct Expr *stop,
           CONST struct Instance *ref,
           struct Instance *rel,
           int *added,
           int i,
           enum relation_errors *err,
           enum find_errors *ferr)
{
  struct relation_term *term = NULL;
  struct gl_list_t *instances;
  unsigned c,len;
  struct Instance *inst;
  struct value_t svalue,cvalue;
  int my_added=0;
  symchar *str;
  CONST struct for_var_t *fvp;  /* for var pointer */
  while (ptr!=stop){
    switch(ExprType(ptr)){
    case e_plus:
    case e_minus:
    case e_times:
    case e_divide:
    case e_power:
    case e_ipower:
    case e_uminus:
      term = CreateOpTerm(ExprType(ptr));
      my_added++;
      AppendTermBuf(term);
      break;
    case e_var:
      str = SimpleNameIdPtr(ExprName(ptr));
      if (str&&TempExists(str)){
        cvalue = TempValue(str);
        switch(ValueKind(cvalue)){
        case integer_value:
          term = CreateIntegerTerm(IntegerValue(cvalue));
          my_added++;
          AppendTermBuf(term);
          break;
        default:
          FPRINTF(ASCERR,"Non-integer temporary variable used in expression.\n");
          *err = incorrect_inst_type;
          term = NULL;
          return 1;
        }
      }else if (GetEvaluationForTable() != NULL && str !=NULL &&
           (fvp=FindForVar(GetEvaluationForTable(),str)) !=NULL ){
        if (GetForKind(fvp)==f_integer){
          term = CreateIntegerTerm(GetForInteger(fvp));
          my_added++;
          AppendTermBuf(term);
        }
        else{
          FPRINTF(ASCERR,
          "Non-integer FOR variable used in expression.\n");
          *err = incorrect_inst_type;
          return 1;
         }
      }
      else{
        instances = FindInstances(ref,ExprName(ptr),ferr);
        if (instances!=NULL){
          if (NextExpr(ptr)==stop){ /* possibly multiple instances */
            len = gl_length(instances);
            for(c=1;c<=len;c++){
              inst = (struct Instance *)gl_fetch(instances,c);
              if ((term=CreateTermFromInst(inst,rel,err))!=NULL){
                AppendTermBuf(term);
                if (my_added++){
                  switch(i){
                  case SUM:
                    term = CreateOpTerm(e_plus);
                    break;
                  case PROD:
                    term = CreateOpTerm(e_times);
                    break;
                  }
                  AppendTermBuf(term);
                }
              }
              else{
                gl_destroy(instances);
                return 1;
              }
            }
            gl_destroy(instances);
          }
          else{         /* single instance */
            if (gl_length(instances)==1){
              inst = (struct Instance *)gl_fetch(instances,1);
              gl_destroy(instances);
              if ((term=CreateTermFromInst(inst,rel,err))!=NULL){
                my_added++;
                AppendTermBuf(term);
              }
              else
                return 1;
              }
              else{
              gl_destroy(instances);
              *err = incorrect_structure;
              FPRINTF(ASCERR,"incorrect_structure in ConvertSubExpr 1\n");
              return 1;
            }
          }
        } else{
          *err = find_error;
          return 1;
        }
      }
      break;
    case e_int:
      term = CreateIntegerTerm(ExprIValue(ptr));
      my_added++;
      AppendTermBuf(term);
      break;
    case e_zero:
      /* this should never happen here */
      term = CreateZeroTerm();
      my_added++;
      AppendTermBuf(term);
      break;
    case e_real:
      term = CreateRealTerm(ExprRValue(ptr),ExprRDimensions(ptr));
      my_added++;
      AppendTermBuf(term);
      break;
    case e_card:
      assert(GetEvaluationContext() == NULL);
      SetEvaluationContext(ref);
      svalue = EvaluateSet(ExprBuiltinSet(ptr),InstanceEvaluateName);
      SetEvaluationContext(NULL);
      cvalue = CardValues(svalue);
      DestroyValue(&svalue);
      switch(ValueKind(cvalue)){
      case integer_value:
        term = CreateIntegerTerm(IntegerValue(cvalue));
        my_added++;
        AppendTermBuf(term);
        break;
      case error_value:
        FigureOutError(cvalue,err,ferr);