/* $Header: /afs/athena.mit.edu/astaff/project/atdev/src/fmax/RCS/functions.c,v 1.7 91/03/05 09:13:27 dot Exp $ */

/*******************************************************************
  Copyright (C) 1990 by the Massachusetts Institute of Technology

   Export of this software from the United States of America is assumed
   to require a specific license from the United States Government.
   It is the responsibility of any person or organization contemplating
   export to obtain such a license before exporting.

WITHIN THAT CONSTRAINT, permission to use, copy, modify, and
distribute this software and its documentation for any purpose and
without fee is hereby granted, provided that the above copyright
notice appear in all copies and that both that copyright notice and
this permission notice appear in supporting documentation, and that
the name of M.I.T. not be used in advertising or publicity pertaining
to distribution of the software without specific, written prior
permission.  M.I.T. makes no representations about the suitability of
this software for any purpose.  It is provided "as is" without express
or implied warranty.

***************************************************************** */

/* Credit to the authors of Gnuplot for some of the functions below.
 * They are all greatly modified from Gnuplot.  See the Gnuplot credit
 * elsewhere in this source
 */

#include <math.h>
#include <stdio.h>
#include <errno.h>
#include "datatypes.h"

extern BOOLEAN undefined;

char *strcpy();

/* prototypes */
Value *pop(), *makecmplx(), *makeint();
double magnitude(), angle(), realpart();

double infnan(iarg)
int iarg;
{
    switch(iarg) {
    case ERANGE:  runerr("Floating point overflow.");
    case -ERANGE: runerr("Division by zero.");
    case EDOM:    runerr("Undefined result.");
    }
}


cnvtointeger(x)
     Value *x;
{
  switch (x->type) {
  case integer:
    return;
  case real:
    x->i.i = (int)x->r.r;
    break;
  case complex:
    x->i.i = (int)x->c.c.r;
    break;
  default:
    internalerr("impossible type conversion");
  }
  x->type = integer;
}

cnvtoreal(x)
     Value *x;
{
  switch (x->type) {
  case integer:
    x->r.r = (double)x->i.i;
    break;
  case real:
    return;
  case complex:
    x->r.r = x->c.c.r;
    break;
  default:
    internalerr("impossible type conversion");
  }
  x->type = real;
}

cnvtocomplex(x)
     Value *x;
{
  switch(x->type) {
  case integer:
    x->c.c.r = (double)x->i.i;
    x->c.c.i =  0.0;
    break;
  case real:
    x->c.c.r = x->r.r;
    x->c.c.i = 0.0;
    break;
  case complex:
    return;
  default:
    internalerr("impossible type conversion");
  }
  x->type = complex;
}

promote(a,b)
     Value *a,*b;
{
  if (a->type == b->type) return;
  
  if (a->type == complex) cnvtocomplex(b);
  else if (b->type == complex) cnvtocomplex(a);
  else if (a->type == real) cnvtoreal(b);
  else if (b->type == real) cnvtoreal(a);
}


/*
static int_check(v)
     Value *v;
{
  if (v->type != integer)
    runerr("non-integer passed to boolean operator");
}
*/
void f_push(),f_pushc(),f_pushd(),f_call(),f_index(),
    f_lnot(),f_bnot(),f_uminus(), f_uplus(),
    f_lor(),f_land(),f_bor(),f_xor(),f_band(),f_eq(),f_ne(),f_gt(),f_lt(),
    f_ge(),f_le(),f_plus(),f_minus(),f_mult(),f_div(),f_mod(),f_power(),
    f_factorial(),f_bool(),f_jump(),f_jumpz(),f_jumpnz(),f_jtern();

void f_real(),f_imag(),f_ang(),f_conjg(),f_sin(),f_cos(),f_tan(),f_asin(),  
       f_acos(),f_atan(),f_sinh(),f_cosh(),f_tanh(),f_int(),f_abs(),f_sgn(),
    f_sqrt(),f_exp(),f_log10(),f_log(),f_besj0(),f_besj1(),f_besy0(),f_besy1(),
    f_gamma(),
    f_floor(),f_ceil(),f_plusover(),f_minusover(),f_multover(), f_divover(),
    f_modover(), f_powerover(), f_lorover(), f_borover(), f_xorover(),
    f_landover(),f_bandover(),f_sizeof(), f_concat();


#define ERR_VAR "undefined variable: "


char *nameofdatatype(d)
datatype d;
{
  switch (d) {
    case integer:   return "integer";
    case real:      return "real"; 
    case complex:   return "complex";
    case array:     return "array";
    case deffunc:   return "user defined function";
    case primfunc:  return "primitive function";
    case undef:     return "undef"; 
    default:        internalerr("nameofdatatype: Unknown type");
  }  
}

static char TMString[80];

TypeMismatch(passed,prmnum,func,expected)
datatype passed;
int prmnum;
char *func;
char *expected;
{
  sprintf(TMString,
	  "Type Mismatch: %s passed as parameter %d to %s; Expecting %s",
	  nameofdatatype(passed),prmnum,func,expected);
  runerr(TMString);
}

TypeMismatchUO(passed,func,expected) /* like above for Unary Operators */
datatype passed;
char *func;
char *expected;
{
  sprintf(TMString,
	  "Type Mismatch: %s applied to %s; Expecting %s",
	  func, nameofdatatype(passed),expected);
  runerr(TMString);
}

TypeMismatchBO(lhs,rhs,func,expected) /* for Binary Operators */
datatype lhs,rhs;
char *func;
char *expected;
{
  sprintf(TMString,
	  "Type Mismatch: %s %s %s; Expecting %s",
	  nameofdatatype(lhs),func,nameofdatatype(rhs),expected);
  runerr(TMString);
}

/* 
 * the two functions below must malloc a chunk of memory, but before they 
 * can use it, they may be blasted back to the command line by errors 
 * when evaluating functions.  Therefore, the pointer to malloc'ed memory
 * is stored in a static variable, and each time we enter the routine,
 * we free any memory we lost last time.  Also, each time we exit normally, 
 * we set the static variable to NULL.  See also NewFuncPlot in plot.c 
 */
/* The above doesn't work, because static variables cannot be used
 * recursively.  I have taken the static stuff out, and this currently
 * constitutes a MEMORY LEAK.  A general scheme for fixing it (also
 * useful in, for example the CalcPoints routine for plotting, is to
 * have a TmpMalloc() function which allocates memory linked in a special
 * list.  After any longjmp() back up to the top, this memory can be freed.
 * There should also be a MakePerm() function to unlink the memory from the
 * temporary list, when any dangerous sections are passed over, and it is
 * guarenteed that no longjmps will be taken.
 */

CalcArray(a,b,f)  /* apply f() to each element of array a and return in b */
Value *a,*b;      /* the caller is responsible for remembering a->vals */
FUNC_PTR f;       /* and free()ing it if desired */
{
  Value *newvals;
  int i;

  if (!isarray(a)) error("CalcArray: not passed array");
  newvals = (Value *) malloc(a->a.N * sizeof(Value));

  for(i=0; i<a->a.N; i++) {          /* we may not exit from this loop */
    f(&(a->a.vals[i]),&(newvals[i]));
  }
  b->type = array;  /* but if we do exit, we are safe now */
  b->a.N = a->a.N;
  b->a.vals = newvals;
}

CalcArrayBO(a,b,c,f)  /* apply f() to pairs of elements in arrays a and b */
Value *a,*b,*c;       /* and return in array c.  See above also */
FUNC_PTR f;
{
  Value *newvals;
  int i;
  Value x;

  if (isarray(a) && isarray(b)) {
    if (a->a.N != b->a.N) 
      runerr("Type Mismatch: arrays have different lengths");
    newvals = (Value *) malloc(a->a.N * sizeof(Value));
    for(i=0; i<a->a.N; i++) {
      f(&(a->a.vals[i]),&(b->a.vals[i]),&(newvals[i]));
    }
    c->type = array;
    c->a.N = a->a.N;
    c->a.vals = newvals;
  }
  else if (isarray(a)) {  /* and b is not an array */
    newvals = (Value *) malloc(a->a.N * sizeof(Value));
    for(i=0; i<a->a.N; i++) {
      x = *b;        /* copy b each time, since f() is allowed to trash args */
      f(&(a->a.vals[i]),&x,&(newvals[i]));
    }
    c->type = array;
    c->a.N = a->a.N;
    c->a.vals = newvals;
  }
  else if (isarray(b))  { /* and a is not an array */
    newvals = (Value *) malloc(b->a.N * sizeof(Value));
    for(i=0; i<b->a.N; i++) {
      x = *a;                     /* copy of a for the same reason */
      f(&x,&(b->a.vals[i]),&(newvals[i]));
    }
    c->type = array;
    c->a.N = b->a.N;
    c->a.vals = newvals;
  }
  else error("CalcArrayBO: neither parameter is an array");
}

CalcOver(a,b,f)
Value *a,*b;
FUNC_PTR f;
{
  Value total;
  int i;

  total = a->a.vals[0];
  for(i=1; i<a->a.N; i++) {
    f(&total,&(a->a.vals[i]),&total);
  }
  *b = total;
}


/* below are a bunch of functions which fmax supports.  They are called
 * with one or two or more Value pointers as arguments and one Value 
 * pointer in which to return the return value.  The address of these 
 * functions will be stored in a table, and they will be called via table
 * lookup.  Since they are part of a rpn system, the value pointers will
 * usually be to the top of the stack.  Often, the return value pointer
 * will be the same as one of the argument value pointers, and the function
 * code must accomodate this.   The functions are allowed to change the
 * values of the passed parameters, so the caller should pass a copy if it
 * wants to preserve the value.  (The functions must modify the parameters
 * before it assigns the return value, since the two pointers may be the 
 * same.)  Currently the only functions which modify their arguments are
 * some of the binary operators which call promote().
 */

void f_lnot(a,b)
     Value *a,*b;
{
  if (a->type == integer) {
    b->type = integer;
    b->i.i = !a->i.i;
  }
  else if (a->type == array) CalcArray(a,b,f_lnot);
  else TypeMismatchUO(a->type,"!","integer");
}


void f_bnot(a,b)
     Value *a,*b;
{
  if (a->type == integer) {
    b->type = integer;
    b->i.i = ~a->i.i;
  }
  else if (a->type == array) CalcArray(a,b,f_bnot);
  else TypeMismatchUO(a->type,"~","integer");
}


void f_bool(a,b)      /* there is no operator associated with this function */
     Value *a,*b;     /* it is just used in the short circuiting for || and && */
{	     
  if (a->type == integer) {
    b->type = integer;
    b->i.i = !!a->i.i;
  }
  else if (a->type == array) CalcArray(a,b,f_bool);
  else TypeMismatchUO(a->type,"???","integer");
}


void f_uminus(a,b)
     Value *a,*b;
{
  switch(a->type) {
  case integer:
    b->type = integer;
    b->i.i = -a->i.i;
    break;
  case real:
    b->type = real;
    b->r.r = -a->r.r;
    break;
  case complex:
    b->type = complex;
    b->c.c.r = -a->c.c.r;
    b->c.c.i = -a->c.c.i;
    break;
  case array:
    CalcArray(a,b,f_uminus);
    break;
  default:
    TypeMismatchUO(a->type,"-","numeric");
  }
}

void f_uplus(a,b)
Value *a,*b;
{
    *b = *a;
}

void f_lor(a,b,c)
     Value *a,*b,*c;
{
  if ((a->type == integer) && (b->type == integer)) {
    c->type = integer;
    c->i.i = a->i.i || b->i.i;
  }
  else if ((a->type == array) || (b->type == array)) CalcArrayBO(a,b,c,f_lor);
  else TypeMismatchBO(a->type,b->type,"||","integers");
}

void f_land(a,b,c)
     Value *a,*b,*c;
{
  if ((a->type == integer) && (b->type == integer)) {
    c->type = integer;
    c->i.i = a->i.i && b->i.i;
  }
  else if ((a->type == array) || (b->type == array)) CalcArrayBO(a,b,c,f_land);
  else TypeMismatchBO(a->type,b->type,"&&","integers");
}


void f_bor(a,b,c)
     Value *a,*b,*c;
{
  if ((a->type == integer) && (b->type == integer)) {
    c->type = integer;
    c->i.i = a->i.i | b->i.i;
  }
  else if ((a->type == array) || (b->type == array)) CalcArrayBO(a,b,c,f_bor);
  else TypeMismatchBO(a->type,b->type,"|","integers");
}


void f_xor(a,b,c)
     Value *a,*b,*c;
{
  if ((a->type == integer) && (b->type == integer)) {
    c->type = integer;
    c->i.i = a->i.i ^ b->i.i;
  }
  else if ((a->type == array) || (b->type == array)) CalcArrayBO(a,b,c,f_xor);
  else TypeMismatchBO(a->type,b->type,"^","integers");
}


void f_band(a,b,c)
     Value *a,*b,*c;
{
  if ((a->type == integer) && (b->type == integer)) {
    c->type = integer;
    c->i.i = a->i.i & b->i.i;
  }
  else if ((a->type == array) || (b->type == array)) CalcArrayBO(a,b,c,f_band);
  else TypeMismatchBO(a->type,b->type,"&","integers");
}



void f_eq(a,b,c) /* note: floating equality is rare because of roundoff error! */
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"==","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_eq);
    return;
  }

  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->i.i = (a->i.i == b->i.i);
    break;
  case real:
    c->i.i = (a->r.r == b->r.r);
    break;
  case complex:
    c->i.i = ((a->c.c.r == b->c.c.r) && (a->c.c.i == b->c.c.i));
  }
  c->type = integer;
  return;
}


void f_ne(a,b,c)
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"!=","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_ne);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->i.i = (a->i.i != b->i.i);
    break;
  case real:
    c->i.i = (a->r.r != b->r.r);
    break;
  case complex:
    c->i.i = ((a->c.c.r != b->c.c.r) || (a->c.c.i != b->c.c.i));
  }
  c->type = integer;
  return;
}

/****************************************************************
 * the following 4 routines order complex numbers by their real
 * components.  This is not necessarily the right thing.  Perhaps
 * we should generate an error when trying to compare them.
 *****************************************************************/

void f_gt(a,b,c)
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,">","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_gt);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->i.i = (a->i.i > b->i.i);
    break;
  case real:
    c->i.i = (a->r.r > b->r.r);
    break;
  case complex:
    c->i.i = (a->c.c.r > b->c.c.r);
  }
  c->type = integer;
  return;
}


void f_lt(a,b,c)
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"<","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_lt);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->i.i = (a->i.i < b->i.i);
    break;
  case real:
    c->i.i = (a->r.r < b->r.r);
    break;
  case complex:
    c->i.i = (a->c.c.r < b->c.c.r);
  }
  c->type = integer;
  return;
}


void f_ge(a,b,c)
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,">=","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_ge);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->i.i = (a->i.i >= b->i.i);
    break;
  case real:
    c->i.i = (a->r.r >= b->r.r);
    break;
  case complex:
    c->i.i = (a->c.c.r >= b->c.c.r);
  }
  c->type = integer;
  return;
}


void f_le(a,b,c)
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"<=","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_le);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->i.i = (a->i.i <= b->i.i);
    break;
  case real:
    c->i.i = (a->r.r <= b->r.r);
    break;
  case complex:
    c->i.i = (a->c.c.r <= b->c.c.r);
  }
  c->type = integer;
  return;
}

/****************************************************************
 * should integer arithmetic be modulo 2**32, or should integer
 * overflows get promoted to reals?
 * idea: use modulo arithmetic, give user a switch to automatically
 *       promote all integers to reals when lexically scanned.
 ****************************************************************/


void f_plus(a,b,c)
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"+","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_plus);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->type = integer;
    c->i.i = a->i.i + b->i.i;
    break;
  case real:
    c->type = real;
    c->r.r = a->r.r + b->r.r;
    break;
  case complex:
    c->type = complex;
    c->c.c.r = a->c.c.r + b->c.c.r;
    c->c.c.i = a->c.c.i + b->c.c.i;
  }
}




void f_minus(a,b,c)
     Value *a,*b,*c;
{
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"-","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_minus);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->type = integer;
    c->i.i = a->i.i - b->i.i;
    break;
  case real:
    c->type = real;
    c->r.r = a->r.r - b->r.r;
    break;
  case complex:
    c->type = complex;
    c->c.c.r = a->c.c.r - b->c.c.r;
    c->c.c.i = a->c.c.i - b->c.c.i;
  }
}


void f_mult(a,b,c)
     Value *a,*b,*c;
{
  register double x,y;
  
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"*","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_mult);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->type = integer;
    c->i.i = a->i.i * b->i.i;
    break;
  case real:
    c->type = real;
    c->r.r = a->r.r * b->r.r;
    break;
  case complex:
    x = a->c.c.r;  y = b->c.c.r;
    c->type = complex;
    c->c.c.r = x*y - a->c.c.i*b->c.c.i;
    c->c.c.i = x*b->c.c.i + a->c.c.i*y;
  }
}


void f_div(a,b,c)
     Value *a,*b,*c;
{
  register double w,x,y,z,square;
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"/","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_div);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    c->type = integer;
    c->i.i = a->i.i / b->i.i;
    break;
  case real:
    c->type = real;
    c->r.r = a->r.r / b->r.r;
    break;
  case complex:
    w = a->c.c.r; x = a->c.c.i;
    y = b->c.c.r; z = b->c.c.i;
    square = y*y + z*z;
    if (square != 0.0) {
      c->type = complex;
      c->c.c.r = (w*y + x*z)/square;
      c->c.c.i = (x*y - w*z)/square;
    }
    else {
      c->type = undef;
      undefined = TRUE;
    }
  }
}


void f_mod(a,b,c)
     Value *a,*b,*c;
{
  if (!((isintorreal(a) || isarray(a)) && (isintorreal(b) || isarray(b)))) 
    TypeMismatchBO(a->type,b->type,"%","integers or real numbers"); 
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_mod);
    return;
  }
  
  if (a->type != b->type) promote(a,b);
  
  if (a->type == integer) {
    c->type = integer;
    c->i.i = a->i.i % b->i.i;
  }
  else {  /* type must == real */
    c->type = real;
    c->r.r = a->r.r - ((int)(a->r.r/b->r.r))*b->r.r;
  }
}


/*****************************************************************
 * NOTE: -x**y parses as -(x**y).
 *       I think this is correct, but check the precedence
 *****************************************************************/

void f_power(a,b,c)
     Value *a,*b,*c;
{
  register int i,t;
  register double mag,ang;
  if (!isnumorarray(a) || !isnumorarray(b)) 
    TypeMismatchBO(a->type,b->type,"**","numbers");  /* never returns */
  if (isarray(a) || isarray(b)) {
    CalcArrayBO(a,b,c,f_power);
    return;
  }
  
  if (a->type != b->type) promote(a,b);  
  
  switch(a->type) {
  case integer:
    if ((b->i.i >= 0) && (b->i.i <= 32)) {
      for(t=1,i=0;i<b->i.i;i++)
	t *= a->i.i;
      c->type = integer;
      c->i.i = t;
      break;  /* get out of case now, else clause drops through */
    }
    else {
      cnvtoreal(a);
      cnvtoreal(b);
    }        /* drop through to next case here */
  case real:
    if ((a->r.r >= 0.0) || (b->r.r == (int)b->r.r)) {
      c->type = real;
      c->r.r = pow(a->r.r,b->r.r);
      break; /* get out not, so we don't drop through */
    }
    else {
      cnvtocomplex(a);
      cnvtocomplex(b);
    }       /* drop through to next case */
  case complex:
    mag = pow(magnitude(a),b->c.c.r);
    ang = angle(a)*b->c.c.r + b->c.c.i;
    c->type = complex;
    c->c.c.r = mag*cos(ang);
    c->c.c.i = mag*sin(ang);
  }
}

static double FactorialTable[] = {  /* size of table is machine dependent */
  1,	                   1,	                   2,
  6,	                  24,	                 120,
  720,	                5040,	               40320,
  362880,	             3628800,	            39916800,
  4790016e+2,	          62270208e+2,	         871782912e+2,
  1307674368e+3,	      20922789888e+3,	     355687428096e+3,
  6402373705728e+3,	  121645100408832e+3,	 243290200817664e+4,
  5109094217170944e+4,	1.1240007277776077e+21,	2.5852016738884977e+22,
  6.2044840173323944e+23,	1.5511210043330986e+25,	4.0329146112660564e+26,
  1.0888869450418352e+28,	3.0488834461171386e+29,	8.841761993739702e+30,
  2.6525285981219106e+32,	8.2228386541779229e+33,	2.6313083693369353e+35,
  8.6833176188118865e+36
};	

void f_factorial(a,b)   /* the numbers 13 and 33 are machine dependent */
     Value *a,*b;
{
  if (isarray(a)) CalcArray(a,b,f_factorial);
  else if (!isinteger(a)) TypeMismatchUO(a->type,"!","integer");
  else if ((a->i.i >= 0) && (a->i.i <= 12)) {
    b->type = integer;
    b->i.i = (int)FactorialTable[a->i.i]; }
  else if ((a->i.i >=13) && (a->i.i <= 33)) {
    b->type = real;
    b->r.r = FactorialTable[a->i.i];
  }
  else {
    b->type = undef;
    undefined = TRUE;
  }
}


/*****************************************************************
 * 
 *    operators above
 *
 *    primitive functions below
 *
 *****************************************************************/





void f_real(a,b)
     Value *a,*b;
{
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = (double) a->i.i;
      break;
    case real: 
      b->type = real;
      b->r.r = a->r.r;
      break;
    case complex:
      b->type = real;
      b->r.r = a->c.c.r;
      break;
    case array:
      CalcArray(a,b,f_real);
      break;
    default:
      TypeMismatch(a->type,1,"real()","numeric");
  }
  return;
}

void f_imag(a,b)
     Value *a,*b;
{
  switch(a->type) {
    case integer:  
    case real: 
      b->type = real;
      b->r.r = 0.0;
      break;
    case complex:
      b->type = real;
      b->r.r = a->c.c.i;
      break;
    case array:
      CalcArray(a,b,f_imag);
      break;
    default:
      TypeMismatch(a->type,1,"imag()","numeric");
  }
  return;
}


void f_ang(a,b)
Value *a,*b;
{
  switch(a->type) {
    case integer: 
      b->type = real;
      b->r.r = (a->i.i > 0)?(0.0):(Pi);
      break;
    case real: 
      b->type = real;
      b->r.r = (a->r.r > 0.0)?(0.0):(Pi);
      break;
    case complex:
      b->type = real;
      b->r.r = angle(a);
      break;
    case array:
      CalcArray(a,b,f_ang);
      break;
    default:
      TypeMismatch(a->type,1,"ang()","numeric");
  }
  return;
}

void f_conjg(a,b)
Value *a,*b;
{
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = (double)a->i.i;
      break;
    case real: 
      b->type = real;
      b->r.r = a->r.r;
      break;
    case complex:
      b->type = complex;
      b->c.c.r = a->c.c.r;
      b->c.c.r = -a->c.c.i;
      break;
    case array:
      CalcArray(a,b,f_conjg);
      break;
    default:
      TypeMismatch(a->type,1,"conjg()","numeric");
  }
  return;
}

void f_sin(a,b)
Value *a,*b;
{
  register double x,y;
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = sin((double)a->i.i);
      break;
    case real: 
      b->type = real;
      b->r.r = sin(a->r.r);
      break;
    case complex:
      b->type = complex;
      x = a->c.c.r;  y = a->c.c.i;
      b->c.c.r = sin(x) * cosh(y);
      b->c.c.i = cos(x) * sinh(y);
      break;
    case array:
      CalcArray(a,b,f_sin);
      break;
    default:
      TypeMismatch(a->type,1,"sin()","numeric");
  }
  return;
}

void f_cos(a,b)
Value *a,*b;
{
  register double x,y;
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = cos((double)a->i.i);
      break;
    case real: 
      b->type = real;
      b->r.r = cos(a->r.r);
      break;
    case complex:
      b->type = complex;
      x = a->c.c.r;  y = a->c.c.i;
      b->c.c.r = cos(x) * cosh(y);
      b->c.c.i = -sin(x) * sinh(y);
      break;
    case array:
      CalcArray(a,b,f_cos);
      break;
    default:
      TypeMismatch(a->type,1,"cos()","numeric");
  }
  return;
}

void f_tan(a,b)
Value *a,*b;
{
  register double den;
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = tan((double)a->i.i);
      break;
    case real: 
      b->type = real;
      b->r.r = tan(a->r.r);
      break;
    case complex:
      den = cos(2*a->c.c.r) + cosh(2*a->c.c.i);
      if (den == 0.0) {
	  b->type = undef;    
	  undefined = TRUE;
      } 
      else {
	  b->type = complex;
	  b->c.c.r = sin(2*a->c.c.r)/den;
	  b->c.c.i = sinh(2*a->c.c.i)/den;
      }
      break;
    case array:
      CalcArray(a,b,f_tan);
      break;
    default:
      TypeMismatch(a->type,1,"tan()","numeric");
  }
  return;
}

void f_asin(a,b)
Value *a,*b;
{
register double alpha, beta, x, y, u, v;
  switch(a->type) {
    case integer:  
    case real: 
      if (a->type == integer) x = (double)a->i.i; else x = a->r.r;
      if (fabs(x) > 1.0) {
	  b->type = undef;    
	  undefined = TRUE;
      }
      else {
	  b->type = real;
	  b->r.r = asin(x);
      }
      break;
    case complex:
      b->type = complex;
      x = a->c.c.r;  y = a->c.c.i;
      u = sqrt((x+1)*(x+1) + y*y); 
      v = sqrt((x-1)*(x-1) + y*y);
      alpha = (u+v)/2;
      beta = (u-v)/2;
      b->c.c.r = asin(beta);
      b->c.c.i = log(alpha + sqrt(alpha*alpha-1));
      break;
    case array:
      CalcArray(a,b,f_asin);
      break;
    default:
      TypeMismatch(a->type,1,"asin()","numeric");
  }
  return;
}

void f_acos(a,b)
Value *a,*b;
{
  register double alpha, beta, x,y,u,v;
  
  switch(a->type) {
    case integer:  
    case real: 
      if (a->type == integer) x = (double)a->i.i; else x = a->r.r;
      if (fabs(x) > 1.0) {
	  b->type = undef;    
	  undefined = TRUE;
      }
      else {
	  b->type = real;
	  b->r.r = acos(x);
      }
      break;
    case complex:
      b->type = complex;
      x = a->c.c.r;  y = a->c.c.i;
      u = sqrt((x+1)*(x+1) + y*y); 
      v = sqrt((x-1)*(x-1) + y*y);
      alpha = (u+v)/2;
      beta = (u-v)/2;
      b->c.c.r = acos(beta);
      b->c.c.i = log(alpha + sqrt(alpha*alpha-1));
      break;
    case array:
      CalcArray(a,b,f_acos);
      break;
    default:
      TypeMismatch(a->type,1,"acos()","numeric");
  }
  return;
}

void f_atan(a,b)
Value *a,*b;
{
register double x, y;
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = atan((double)a->i.i);
      break;
    case real: 
      b->type = real;
      b->r.r = atan(a->r.r);
      break;
    case complex:
      x = a->c.c.r;  y = a->c.c.i;
      if ((x == 0.0) && (fabs(y) == 1.0)) {
	  b->type = undef;    
	  undefined = TRUE;
      } 
      else {
	  b->type = complex;
	  b->c.c.r = atan(2*x/(1-x*x-y*y));
	  b->c.c.i = log((x*x+(y+1)*(y+1))/(x*x+(y-1)*(y-1)))/4;
      }
      break;
    case array:
      CalcArray(a,b,f_atan);
      break;
    default:
      TypeMismatch(a->type,1,"atan()","numeric");
  }
  return;
}

void f_sinh(a,b)
Value *a,*b;
{
  register double x,y;
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = sinh((double)a->i.i);
      break;
    case real: 
      b->type = real;
      b->r.r = sinh(a->r.r);
      break;
    case complex:
      b->type = complex;
      x = a->c.c.r;  y = a->c.c.i;
      b->c.c.r = sinh(x) * cos(y);
      b->c.c.i = cosh(x) * sin(y);
      break;
    case array:
      CalcArray(a,b,f_sinh);
      break;
    default:
      TypeMismatch(a->type,1,"sinh()","numeric");
  }
  return;
}

void f_cosh(a,b)
Value *a,*b;
{
  register double x,y;
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = cosh((double)a->i.i);
      break;
    case real: 
      b->type = real;
      b->r.r = cosh(a->r.r);
      break;
    case complex:
      b->type = complex;
      x = a->c.c.r;  y = a->c.c.i;
      b->c.c.r = cosh(x) * cos(y);
      b->c.c.i = sinh(x) * sin(y);
      break;
    case array:
      CalcArray(a,b,f_sinh);
      break;
    default:
      TypeMismatch(a->type,1,"cosh()","numeric");
  }
  return;
}

void f_tanh(a,b)
Value *a,*b;
{
  register double den;
  switch(a->type) {
    case integer:  
      b->type = real;
      b->r.r = tanh((double)a->i.i);
      break;
    case real: 
      b->type = real;
      b->r.r = tanh(a->r.r);
      break;
    case complex:
      b->type = complex;
      den = cosh(2*a->c.c.r) + cos(2*a->c.c.i);
      b->c.c.r = sinh(2*a->c.c.r)/den;
      b->c.c.i = sin(2*a->c.c.i)/den;
      break;
    case array:
      CalcArray(a,b,f_tanh);
      break;
    default:
      TypeMismatch(a->type,1,"tanh()","numeric");
  }
  return;
}

void f_int(a,b)
Value *a,*b;
{
  switch(a->type) {
    case integer:  
      b->type = integer;
      b->i.i =  a->i.i;
      break;
    case real: 
      b->type = integer;
      b->i.i = (int)a->r.r;
      break;
    case complex:
      b->type = integer;
      b->i.i = (int)a->c.c.r;
      break;
    case array:
      CalcArray(a,b,f_int);
      break;
    default:
      TypeMismatch(a->type,1,"int()","numeric");
  }
  return;
}


void f_abs(a,b)
Value *a,*b;
{
  switch (a->type) {
    case integer:
      b->type = integer;
      b->i.i = abs(a->i.i);
      break;
    case real:
      b->type = real;
      b->r.r = fabs(a->r.r);
      break;
    case complex:
      b->type = real;
      b->r.r = sqrt(a->c.c.r*a->c.c.r + a->c.c.i*a->c.c.i);
      break;
    case array:
      CalcArray(a,b,f_abs);
      break;
    default:
      TypeMismatch(a->type,1,"abs()","numeric");
  }
  return;
}

void f_sgn(a,b)
Value *a,*b;
{
  switch (a->type) {
    case integer:
      b->type = integer;
      b->i.i = (a->i.i > 0) ? 1 : ((a->i.i < 0) ? -1 : 0);
      break;
    case real:
      b->type = integer;
      b->r.r = (a->r.r > 0.0) ? 1 : ((a->r.r < 0.0) ? -1 : 0);
      break;
    case complex:
      b->type = integer;
      b->r.r = (a->c.c.r > 0.0) ? 1 : ((a->c.c.r < 0.0) ? -1 : 0);
      break;
    case array:
      CalcArray(a,b,f_sgn);
      break;
    default:
      TypeMismatch(a->type,1,"sgn()","numeric");
  }
  return;
}


void f_sqrt(a,b)
Value *a,*b;
{
  register double x,mag,ang;
  switch(a->type) {
    case integer:  
    case real: 
      if (a->type == integer) x = (double)a->i.i; else x = a->r.r;
      if (x >= 0.0) {
	  b->type = real;
	  b->r.r = sqrt(x);
      }
      else {
	  b->type = complex;
	  b->c.c.r = 0.0;
	  b->c.c.i = sqrt(-x);
      }
      break;
    case complex:
      b->type = complex;
      mag = sqrt(magnitude(a));
      ang = angle(a);
      if (ang < 0.0) ang += 2*Pi;
      ang /= 2;
      b->c.c.r = mag*cos(ang);
      b->c.c.i = mag*sin(ang);
      break;
    case array:
      CalcArray(a,b,f_sqrt);
      break;
    default:
      TypeMismatch(a->type,1,"sqrt()","numeric");
  }
  return;
}


void f_exp(a,b)
Value *a,*b;
{
  register double mag;
  switch(a->type) {
    case integer:
      b->type = real;
      b->r.r = exp((double)a->i.i);
      break;
    case real:
      b->type = real;
      b->r.r = exp(a->r.r);
      break;
    case complex:
      b->type = complex;
      mag = exp(a->c.c.r);
      b->c.c.r = mag*cos(a->c.c.i);
      b->c.c.i = mag*sin(a->c.c.i);
      break;
    case array:
      CalcArray(a,b,f_exp);
      break;
    default:
      TypeMismatch(a->type,1,"exp()","numeric");
  }
  return;
}


void f_log10(a,b)
Value *a,*b;
{
  register double x,mag,ang;

  switch(a->type) {
    case integer:
    case real:
      if (a->type == integer) x = (double)a->i.i; else x = a->r.r;
      if (x > 0.0) {
	  b->type = real;
	  b->r.r = log(x)/2.3025848865509033;
      }
#ifdef NOTDEF
/* 11/29/90 Dorothy Bowe - Can't take the log of a negative number */
      else {
	  b->type = complex;
	  b->c.c.r = log(x)/2.3025848865509033;
	  b->c.c.i = Pi/2.3025848865509033;
      }
#endif
      break;
    case complex:
      b->type = complex;
      mag = magnitude(a);  ang = angle(a);
      b->c.c.r = log(mag)/2.3025848865509033;
      b->c.c.i = ang/2.3025848865509033;
      break;
    case array:
      CalcArray(a,b,f_log10);
      break;
    default:
      TypeMismatch(a->type,1,"log10()","numeric");
  }
  return;
}


void f_log(a,b)
Value *a,*b;
{
  register double x,mag,ang;
  switch(a->type) {
    case integer:
    case real:
      if (a->type == integer) x = (double)a->i.i; else x = a->r.r;
      if (x > 0.0) {
	  b->type = real;
	  b->r.r = log(x);
      }
      else {
	  b->type = complex;
	  b->c.c.r = log(-x);
	  b->c.c.i = Pi;
      }
      break;
    case complex:
      b->type = complex;
      mag = magnitude(a);  ang = angle(a);
      b->c.c.r = log(mag);
      b->c.c.i = ang;
      break;
    case array:
      CalcArray(a,b,f_log);
      break;
    default:
      TypeMismatch(a->type,1,"log()","numeric");
  }
  return;
}


void f_besj0(a,b)   
Value *a,*b;
{
     float z;
     double xx,y,ans1,ans2,val;
     Value c;
     
     switch (a->type) {
	case integer:
	case real:
	case complex:
	  f_abs (a,&c);
	  val = (a->type == integer) ? a->i.i : a->r.r;
	  if (((c.type == integer) && (c.i.i < 8)) ||
	        ((c.type == real) && (c.r.r < 8.0))) {
	       y = val * val;
	       ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
                        +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
	       ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
			+y*(59272.64853+y*(267.8532712+y*1.0))));
	       b->type = real;
	       b->r.r = ans1/ans2;
	  }
	  else {
	       b->type = real;
	       if (c.type == integer) 
		    b->r.r = c.i.i;
	       else
		    b->r.r = c.r.r;
	       z = 8.0 / b->r.r;
	       y = z * z;
	       xx = b->r.r - 0.785398164;
	       ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
		        +y*(-0.2073370639e-5+y*0.2093887211e-6)));
	       ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			-y*0.934935152e-7)));
	       b->r.r = sqrt(0.636619772/b->r.r)*
			 (cos(xx)*ans1-z*sin(xx)*ans2);
	  }
	  break;
	case array:
	  CalcArray(a,b,f_besj0);
	  break;
	default:
	  TypeMismatch(a->type,1,"besj0()","numeric");
     }
     return;
}

void f_besj1(a,b)
Value *a,*b;
{
     float z;
     double xx,y,ans1,ans2,val;
     Value c;

     switch (a->type) {
	case integer:
	case real:
	case complex:
	  f_abs (a,&c);
	  b->type = real;
	  val = (a->type == integer) ? a->i.i : a->r.r;
	  if (((c.type == integer) && (c.i.i < 8)) ||
	        ((c.type == real) && (c.r.r < 8.0))) {
	       y = val * val;
	       ans1 = val*(72362614232.0+y*(-7895059235.0+y*(242396853.1
			    +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
	       ans2 = 144725228442.0+y*(2300535178.0+y*(18583304.74
			+y*(99447.43394+y*(376.9991397+y*1.0))));
	       b->r.r = ans1/ans2;
	  }
	  else {
	       if (c.type == integer)
		    b->r.r = c.i.i;
	       else
		    b->r.r = c.r.r;
	       z = 8.0 / b->r.r;
	       xx = b->r.r - 2.356194491;
	       y = z * z;
	       ans1 = 1.0+y*(0.183105e-2+y*(-0.3516396496e-4
		      +y*(0.2457520174e-5+y*(-0.240337019e-6))));
	       ans2 = 0.04687499995+y*(-0.2002690873e-3
		      +y*(0.8449199096e-5+y*(-0.88228987e-6
		      +y*0.105787412e-6)));
	       b->r.r = sqrt(0.636619772/b->r.r)*
			 (cos(xx)*ans1-z*sin(xx)*ans2);
	       if (val < 0)
		    b->r.r = -b->r.r;
	       }
	  break;
	case array:
	  CalcArray(a,b,f_besj1);
	  break;
	default:
	  TypeMismatch(a->type,1,"besj1()","numeric");
     }
     return;
}


void f_besy0(a,b)   /* y0(a) = -cos(a)/a */
Value *a,*b;
{
     float z;
     double xx,y,ans1,ans2;
     void f_besj0();
     Value c,x;
     double val;


     switch (a->type) {
	case integer:
	case real:
	case complex:
	  f_abs (a,&c);
	  b->type = real;
	  val = (a->type == integer) ? a->i.i : a->r.r;
	  if (((c.type == integer) && (c.i.i < 8)) ||
	        ((c.type == real) && (c.r.r < 8.0))) {
	       y = val * val;
	       ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
			+y*(10879881.29+y*(-86327.92757+y*228.4622733))));
	       ans2=40076544269.0+y*(745249964.8+y*(7189466.438
			+y*(47447.26470+y*(226.1030244+y*1.0))));
	       f_besj0(a,&x);
	       if (val == 0)
		    val = HUGE;
	       else
		    val = (a->type == integer) ? log(a->i.i) : log(a->r.r);
	       if (x.type == integer)
		    b->r.r = (ans1/ans2)+0.636619772*x.i.i*val;
	       else
		    b->r.r = (ans1/ans2)+0.636619772*x.r.r*val;
	  }
	  else {
	       b->r.r = (c.type == integer) ? c.i.i : c.r.r;
	       if (val == 0)
		    z = 0;
	       else
		    z = 8.0 / val;
	       xx = val - 0.785398164;
	       y=z*z;
	       ans1 = 1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			   +y*(-0.2073370639e-5+y*0.2093887211e-6)));
	       ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			   +y*(-0.6911147651e-5+y*(0.7621095161e-6
			   +y*(-0.934945152e-7))));
	       b->r.r = sqrt(0.636619772/val)*(sin(xx)*ans1+z*cos(xx)*ans2);
	  }
	  break;
	case array:
	  CalcArray(a,b,f_besy0);
	  break;
	default:
	  TypeMismatch(a->type,1,"besy0()","numeric");
     }
     return;

}

void f_besy1(a,b)
Value *a,*b;
{
     float z;
     double xx,y,ans1,ans2,val;
     void f_besj1();
     Value c,x;

     switch (a->type) {
	case integer:
	case real:
	case complex:
	  f_abs (a,&c);
	  b->type = real;
	  val = (a->type == integer) ? a->i.i : a->r.r;
	  if (((c.type == integer) && (c.i.i < 8)) ||
	      ((c.type == real) && (c.r.r < 8.0))) {
		y = val * val;
		ans1 = val*(-0.4900604943e13+y*(0.1275274390e13
			+y*(-0.5153438139e11+y*(0.7349264551e9
			+y*(-0.4237922726e7+y*0.8511937935e4)))));
		ans2 = 0.2499580570e14+y*(0.4244419664e12
			+y*(0.3733650367e10+y*(0.2245904002e8
			+y*(0.1020426050e6+y*(0.3549632885e3+y)))));
		f_besj1(a, &x);
		if (val == 0)
		     val = HUGE;
		if (x.type == integer)
		     b->r.r = (ans1/ans2)+0.636619772*
			  (x.r.r)*log(val)-1.0/val;
		else
		     b->r.r = (ans1/ans2)+0.636619772*
			  (x.r.r)*log(val)-1.0/val;
	   }
	  else {
	       if (val == 0)
		    z = 0;
	       else
		    z = 8.0 / val;
	       xx = val - 2.356194491;
	       y=z*z;
	       ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
			+y*(0.2457520174e-5+y*(-0.240337019e-6))));
	       ans2=0.04687499995+y*(-0.2002690873e-3
			+y*(0.8449199096e-5+y*(-0.88228987e-6
			+y*0.105787412e-6)));
		b->r.r=sqrt(0.636619772/val)*(sin(xx)*ans1+z*cos(xx)*ans2);
	  }
	  break;
	case array:
	  CalcArray(a,b,f_besy1);
	  break;
	default:
	  TypeMismatch(a->type,1,"besy1()","numeric");
     }
	  
}


void f_floor(a,b)
Value *a,*b;
{
  switch (a->type) {
    case integer: 
      b->type = integer;
      b->i.i = a->i.i;
      break;
    case real:
      b->type = integer;
      b->i.i = (int)floor(a->r.r);
      break;
    case complex:
      b->type = complex;
      b->c.c.r = floor(a->c.c.r);
      b->c.c.i = floor(a->c.c.i);
      break;
    case array:
      CalcArray(a,b,f_floor);
      break;
    default:
      TypeMismatch(a->type,1,"floor()","numeric");
  }
  return;
}


void f_ceil(a,b)
Value *a,*b;
{
  switch (a->type) {
    case integer: 
      b->type = integer;
      b->i.i = a->i.i;
      break;
    case real:
      b->type = integer;
      b->i.i = (int)ceil(a->r.r);
      break;
    case complex:
      b->type = complex;
      b->c.c.r = ceil(a->c.c.r);
      b->c.c.i = ceil(a->c.c.i);
      break;
    case array:
      CalcArray(a,b,f_ceil);
      break;
    default:
      TypeMismatch(a->type,1,"ceil()","numeric");
  }
  return;
}


void f_gamma(a,b)
Value *a,*b;
{
  extern int signgam;
  register double y;

  switch(a->type) {
    case integer:
    case real:
      if (a->type == integer) y = (double)a->i.i; else y = a->r.r;
      y = gammln(y);
      if (y > 88.0) {
	b->type = undef;    
	undefined = TRUE;
      }
      else {
	b->type = real;
	b->r.r = signgam * exp(y);
      }
      break;
    case array:
      CalcArray(a,b,f_gamma);
      break;
    default:
      TypeMismatch(a->type,1,"gamma()","integer or real");
      break;
  }
  return;
}


void f_plusover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"+/","array");
  CalcOver(a,b,f_plus);
}

void f_minusover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"-/","array");
  CalcOver(a,b,f_minus);
}

void f_multover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"*/","array");
  CalcOver(a,b,f_mult);
}

void f_divover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"//","array");
  CalcOver(a,b,f_div);
}

void f_modover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"%/","array");
  CalcOver(a,b,f_mod);
}

void f_powerover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"**/","array");
  CalcOver(a,b,f_power);
}

void f_lorover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"||/","array");
  CalcOver(a,b,f_lor);
}

void f_borover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"|/","array");
  CalcOver(a,b,f_bor);
}

void f_xorover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"^/","array");
  CalcOver(a,b,f_xor);
}

void f_landover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"&&/","array");
  CalcOver(a,b,f_land);
}

void f_bandover(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"&/","array");
  CalcOver(a,b,f_band);
}

void f_sizeof(a,b)
Value *a,*b;
{
  if (!isarray(a)) TypeMismatchUO(a->type,"$","array");
  b->type = integer;
  b->i.i = a->a.N;
}

void f_concat(a,b,c)
Value *a,*b,*c;
{
    int count = 0;
    int i;

    if (a->type == array) count += a->a.N; else count += 1;
    if (b->type == array) count += b->a.N; else count += 1;
    c->type = array;
    c->a.N = c->a.maxsize = count;
    c->a.vals = (Value *)malloc(count * sizeof(Value));

    count = 0;
    if (a->type == array) {
	for(i=0; i < a->a.N; i++) c->a.vals[count++] = a->a.vals[i];
    }
    else
	c->a.vals[count++] = *a;

    if (b->type == array) {
	for(i=0; i < b->a.N; i++) c->a.vals[count++] = b->a.vals[i];
    }
    else
	c->a.vals[count++] = *b;
    
}