#if ED				// see autodif.h

static char **_IndepName=NULL;	// List of indep vars
static int _IndepNum=0;

static FILE *EDStream_;		// Output for External Differentiator

#endif	/* ED */

int _Niv_;			// Number of Independent Variables
int _Hod_;			// Highest order of derivatives

#define MAXDEP_		10	// max static depth of calls to u.d.f.
#define MAXCALLS_	256	// max number of static calls to u.d.f.

typedef char Set[MAXCALLS_/8];

typedef struct {
  int Max;
  Set CurSet;
  Set OrSet;
} SetOfIndex, *SetOfIndexPtr;

static int FuncLevel_;		// Level of u.d.f. calls

static SetOfIndex sf_set[MAXDEP_];	// reserved indices to sv__ and sd__
static SetOfIndex uf_set[MAXDEP_];	// reserved indices to value_ and grad_

static void ResetIndex(void);
static int sfGetIndex(void);

#define Der0		Expr
#define GetDer0_	GetExpr_

// A bit of psudo-science.
// Let _Niv_ be the number of independent variables v(i), i=0,...,_Niv_-1.
// The map ind2_(i,j)=i+j*(j+1)/2 maps indices of the variables to the
//       interval [0..ind2(_Niv_-1,_Niv_-1)] provided that i<=j.
// The map ind3_(i,j,k)=i+j*(j+1)/2+k*(k+1)*(k+2)/6 maps indices of the variables to the
//       interval [0..ind3(_Niv_-1,_Niv_-1,_Niv_-1)] provided that i<=j<=k.

// Maps subscripts of indep vars to a subscript of appropriate second derivative
static int ind2(int i, int j) {
  return ind2_(i,j);
} 

// Maps subscripts of indep vars to a subscript of appropriate third derivative
static int ind3(int i, int j, int k) {
  return ind3_(i,j,k);
} 

#if ED

static int ind4(int i, int j, int k, int l) {
  return ind4_(i,j,k,l);
} 

static int ind5(int i, int j, int k, int l, int m) {
  return ind5_(i,j,k,l,m);
} 

#endif	/* ED */

//========================================================
// Constructors

// NB: varno can be nonzero only if kind==INDEP
Gradient::Gradient(char *name, int varno, int vardim, VarKind kind) {
  char ***p;
  int i;
  char buf[100];
  Index=0;
  Vardim=vardim;
  Kind=kind;
  Parent=(Gradient *)(name ? this : NULL);
  Expr=(char **)MemNewN((int)(vardim*sizeof(char *)));
  *buf='\0';
  if (vardim>1)
    for (i=0; i<vardim; i++) {
      if (name) sprintf(buf,"%s[%i]",name,(int)i);
      Expr[i]=(char *)StringSave(buf);
    }
  else Expr[0]=(char *)StringSave(name);
#if ED
  if (kind==INDEP && _Hod_>3) {
    /* Create list of names of independent variables */
    char **p;
    int i;
    if (_IndepNum<varno+vardim-1) {
      if (_IndepNum) {
	p=(char**)MemNew((varno+vardim-1)*sizeof(char *));
	MemCopy(p,_IndepName,_IndepNum*sizeof(char *));
	for (i=0; i<_IndepNum; i++) p[i]=_IndepName[i];
	MemFree(_IndepName);
	_IndepName=p;
      } else {
	_IndepName=(char**)MemNew((varno+vardim-1)*sizeof(char *));
      }
    _IndepNum=varno+vardim-1;
    }
    for (i=0; i<vardim; i++) {
      MemFree(_IndepName[varno-1+i]);
      _IndepName[varno-1+i]=(char*)StrSave(Expr[i]);
    }
  }
#endif	/* ED */
  switch (_Hod_) {
    case 3:
      n3=ind3(_Niv_-1,_Niv_-1,_Niv_-1)+1;
      Der3=(char ***)MemNew((int)((n3+1)*vardim*sizeof(char *)));
      for (p=Der3+vardim, i=0; i<vardim; p+=n3, i++)
	Der3[i]=(char **)p;
    case 2:
      n2=ind2(_Niv_-1,_Niv_-1)+1;
      Der2=(char ***)MemNew((int)((n2+1)*vardim*sizeof(char *)));
      for (p=Der2+vardim, i=0; i<vardim; p+=n2, i++)
	Der2[i]=(char **)p;
    case 1:
      Der1=(char ***)MemNew((int)((_Niv_+1)*vardim*sizeof(char *)));
      for (p=Der1+vardim, i=0; i<vardim; p+=_Niv_, i++) {
	Der1[i]=(char **)p;
	if (varno)
	  Der1[i][varno-1+i]=(char *)StrSave("1");	/* others are yet NULLs */
      }
      break;
#if ED
    case 4:
    case 5:
      // Maple: check for names reserved by Maple
      if (name && sscanf(name,"t%i",&i)==1)
	fprintf(MStream_,"#error Name '%s' is reserved by Maple\n",name);
      return;
    default:
      fprintf(MStream_,"#error *** _Hod_>5 || _Hod_<1 ***\n");
#else	/* ED */
    default:
      fprintf(MStream_,"#error *** _Hod_>3 || _Hod_<1 ***\n");
#endif	/* ED */
  }
  if (kind==TEMP) {
    fprintf(MStream_,"  double");
    switch (_Hod_) {
      case 3:
	if (vardim>1)
	  fprintf(MStream_," d3_%s[%i][%i],",name,vardim,n3);
	else
	  fprintf(MStream_," d3_%s[%i],",name,n3);
      case 2:
	if (vardim>1)
	  fprintf(MStream_," d2_%s[%i][%i],",name,vardim,n2);
	else
	  fprintf(MStream_," d2_%s[%i],",name,n2);
      case 1:
	if (vardim>1)
	  fprintf(MStream_," %s[%i],d1_%s[%i][%i];\n",name,vardim,name,vardim,_Niv_);
	else
	  fprintf(MStream_," %s,d1_%s[%i];\n",name,name,_Niv_);
    }
  }
}

// f is a double const
Gradient::Gradient(double c) {
  char buf[20];
  Index=0;
  Vardim=1;
  Kind=DEF;
  Parent=NULL;
  Expr=(char **)MemNewN(sizeof(char *));
  sprintf(buf,"%g",c);
  Expr[0]=(char *)StringSave(buf);
  switch (_Hod_) {
    case 3:
      n3=ind3(_Niv_-1,_Niv_-1,_Niv_-1)+1;
      Der3=(char ***)MemNew((int)((n3+1)*sizeof(char *)));
      Der3[0]=(char **)Der3+1;
    case 2:
      n2=ind2(_Niv_-1,_Niv_-1)+1;
      Der2=(char ***)MemNew((int)((n2+1)*sizeof(char *)));
      Der2[0]=(char **)Der2+1;
    case 1:
      Der1=(char ***)MemNew((int)((_Niv_+1)*sizeof(char *)));
      Der1[0]=(char **)Der1+1;
      break;
    case 4:
    case 5:
      break;
#if ED
    default:
      fprintf(MStream_,"#error *** _Hod_>5 || _Hod_<1 ***\n");
#else	/* ED */
    default:
      fprintf(MStream_,"#error *** _Hod_>3 || _Hod_<1 ***\n");
#endif	/* ED */
  }
}

// f is an int const
Gradient::Gradient(int i) {
  char buf[20];
  Index=0;
  Vardim=1;
  Kind=DEF;
  Parent=NULL;
  Expr=(char **)MemNewN(sizeof(char *));
  sprintf(buf,"%i",i);
  Expr[0]=(char *)StringSave(buf);
  switch (_Hod_) {
    case 3:
      n3=ind3(_Niv_-1,_Niv_-1,_Niv_-1)+1;
      Der3=(char ***)MemNew((int)((n3+1)*sizeof(char *)));
      Der3[0]=(char **)Der3+1;
    case 2:
      n2=ind2(_Niv_-1,_Niv_-1)+1;
      Der2=(char ***)MemNew((int)((n2+1)*sizeof(char *)));
      Der2[0]=(char **)Der2+1;
    case 1:
      Der1=(char ***)MemNew((int)((_Niv_+1)*sizeof(char *)));
      Der1[0]=(char **)Der1+1;
      break;
    case 4:
    case 5:
      break;
#if ED
    default:
      fprintf(MStream_,"#error *** _Hod_>5 || _Hod_<1 ***\n");
#else	/* ED */
    default:
      fprintf(MStream_,"#error *** _Hod_>3 || _Hod_<1 ***\n");
#endif	/* ED */
  }
}

// Copy constructor
Gradient::Gradient(const Gradient &g) {
  char ***p;
  int i,j;
  Vardim=g.Vardim;
  Index=g.Index;
  Kind=g.Kind;
  Parent=g.Parent;
  n2=g.n2;
  n3=g.n3;
  if (g.Expr) {
    Expr=(char **)MemNewN((int)(g.Vardim*sizeof(char *)));
    for (i=0; i<g.Vardim; i++)
      Expr[i]=(char *)StringSave(g.Expr[i]);
  } else Expr=NULL;
  if (g.Der1) {
    switch (_Hod_) {
      case 3:
	Der3=(char ***)MemNewN((int)((n3+1)*g.Vardim*sizeof(char *)));
	for (p=Der3+g.Vardim, i=0; i<g.Vardim; p+=n3, i++) {
	  Der3[i]=(char **)p;
	  for (j=0; j<n3; j++)
	    ((char**)p)[j]=(char *)StringSave(g.Der3[i][j]);
	}
      case 2:
	Der2=(char ***)MemNewN((int)((n2+1)*g.Vardim*sizeof(char *)));
	for (p=Der2+g.Vardim, i=0; i<g.Vardim; p+=n2, i++) {
	  Der2[i]=(char **)p;
	  for (j=0; j<n2; j++)
	    ((char**)p)[j]=(char *)StringSave(g.Der2[i][j]);
	}
      case 1:
	Der1=(char ***)MemNewN((int)((_Niv_+1)*g.Vardim*sizeof(char *)));
	for (p=Der1+g.Vardim, i=0; i<g.Vardim; p+=_Niv_, i++) {
	  Der1[i]=(char **)p;
	  for (j=0; j<_Niv_; j++)
	    ((char**)p)[j]=(char *)StringSave(g.Der1[i][j]);
	}
    }
  } else Der1=NULL;
}

// Destructor
Gradient::~Gradient() {
  int i,j;
  if (Expr) {
    for (i=0; i<Vardim; i++) MemFree(Expr[i]);
    MemFree(Expr);
  }
  if (Der1) {
    switch (_Hod_) {
      case 3:
	for (i=0; i<Vardim; i++) {
	  for (j=0; j<n3; j++) MemFree(Der3[i][j]);
	}
	MemFree(Der3);
      case 2:
	for (i=0; i<Vardim; i++) {
	  for (j=0; j<n2; j++) MemFree(Der2[i][j]);
	}
	MemFree(Der2);
      case 1:
	for (i=0; i<Vardim; i++) {
	  for (j=0; j<_Niv_; j++) MemFree(Der1[i][j]);
	}
	MemFree(Der1);
    }
  }
}

static char *Mult(char *p, char *q, char freeptrs);


//====== A d d ==================================

static char *Add(char *p, char *q, char freeptrs) {
  char *r;
  int len;
  if (p) {
    if (q) {
      len=(int)StrLen(p)+(int)StrLen(q)+1;
      if (*q=='-') {
	r=(char *)MemNewN(len);
	sprintf(r,"%s%s",p,q);
      } else
	if (StrCmp(p,q)) {
	  r=(char *)MemNewN(len+1);
	  sprintf(r,"%s+%s",p,q);
	} else {
	  r=Mult("2",p,0);
	}
    } else {
      r=(char *)StringSave(p);
    }
  } else {
    if (q) {
      r=(char *)StringSave(q);
    } else {
      r=NULL;
    }
  }
  switch (freeptrs) {
    case 'l':
      MemFree(p);
      break;
    case 'r':
      MemFree(q);
      break;
    case 'b':
      MemFree(p);
      MemFree(q);
  }
  return r;
}

// f+g
Gradient operator+(const Gradient &f, const Gradient &g) {
  int j;
  Gradient t;
  // construct expression f+g
  t.Expr[0]=Add(f.Expr[f.Index],g.Expr[g.Index],0);
  switch (_Hod_) {
    case 3:
      // construct (f+g)'''=f'''+g'''
      for (j=0; j<f.n3; j++)
	t.Der3[0][j]=Add(f.Der3[f.Index][j],g.Der3[g.Index][j],0);
    case 2:
      // construct (f+g)''=f''+g''
      for (j=0; j<f.n2; j++)
	t.Der2[0][j]=Add(f.Der2[f.Index][j],g.Der2[g.Index][j],0);
    case 1:
      // construct (f+g)'=f'+g'
      for (j=0; j<_Niv_; j++)
	t.Der1[0][j]=Add(f.Der1[f.Index][j],g.Der1[g.Index][j],0);
  }
  return t;
}

// f+const
Gradient operator+(const Gradient &f, double c) {
  Gradient t(c);
  return f+t;
}

// const+f
Gradient operator+(double c, const Gradient &g) {
  Gradient t(c);
  return t+g;
}

// f+intconst
Gradient operator+(const Gradient &f, int i) {
  Gradient t((double)i);
  return f+t;
}

// intconst+f
Gradient operator+(int i, const Gradient &g) {
  Gradient t((double)i);
  return t+g;
}

//====== M u l t i p l y ========================

static char EqualToOne(char *p) {
  if (*p++!='1') return 0;
  switch (*p++) {
    case 0: return 1;
    case '.': return *p++=='0' ? !*p : 0;
    default: return 0;
  }
}

#define PRTY_UNARY	0
#define PRTY_MULT	1
#define PRTY_ADD	2
#define PRTY_REL	3
#define PRTY_EQU	4
#define PRTY_AND	5
#define PRTY_OR		6
#define PRTY_COND	7
#define PRTY_ASS	8
#define PRTY_		PRTY_ASS+1

static char *PrtyUnary[]={"!","-","+",0};
static char lPrtyUnary[]={  1,  1,  1};
static char *PrtyMult[]={"*","/",0};
static char lPrtyMult[]={  1,  1};
static char *PrtyAdd[]={"+","-",0};
static char lPrtyAdd[]={  1,  1};
static char *PrtyRel[]={">","<",">=","<=",0};
static char lPrtyRel[]={  1,  1,   2,   2};
static char *PrtyEqu[]={"==","!=",0};
static char lPrtyEqu[]={   2,   2};
static char *PrtyAnd[]={"&&",0};
static char lPrtyAnd[]={   2};
static char *PrtyOr[]={"||",0};
static char lPrtyOr[]={   2};
static char *PrtyCont[]={"?",0};
static char lPrtyCont[]={  1};
static char *PrtyAss[]={"=","+=","-=","*=","/=",0};
static char lPrtyAss[]={  1,   2,   2,   2,   2};
static char **Prty[PRTY_]={
  PrtyUnary,
  PrtyMult,
  PrtyAdd,
  PrtyRel,
  PrtyEqu,
  PrtyAnd,
  PrtyOr,
  PrtyCont,
  PrtyAss
};
static char *lPrty[PRTY_]={
  lPrtyUnary,
  lPrtyMult,
  lPrtyAdd,
  lPrtyRel,
  lPrtyEqu,
  lPrtyAnd,
  lPrtyOr,
  lPrtyCont,
  lPrtyAss
};


static int NeedParentheses(char *p, int priority) {
  int i,len,level,prty,o;
  char c;
  for (i=level=0, len=(int)StringLen(p); i<len; i++) {
    if (isalnum(c=p[i])) continue; 
    switch (c) {
      case '[':
      case '(': level++; break;
      case ']':
      case ')': level--; break;
      default:
	if (level==0)
	  for (prty=priority; prty<PRTY_; prty++)
	    for (o=0; Prty[prty][o]; o++) {
	      switch (lPrty[prty][o]) {
		case 1:
		  if (c==Prty[prty][o][0]) return 2;
		  break;
		case 2:
		  if (c==Prty[prty][o][0] && p[i+1]==Prty[prty][o][1]) return 2;
		  break;
	      }
	    }
    }
  }
  return 0;
}

static int IndeX[3];	// indices of independent variables (i<=j<=k)
#define i IndeX[0]
#define j IndeX[1]
#define k IndeX[2]

static char *Interpret(const Gradient &f, const Gradient &g);
static char *desc;

static char *Mult(char *p, char *q, char freeptrs) {
  char *r;
  char *f;
  int lneed,rneed;
  if (p && q) {
    if (EqualToOne(p)) r=(char *)StringSave(q);
    else {
      if (EqualToOne(q)) r=(char *)StringSave(p);
      else {
	lneed=NeedParentheses(p,PRTY_ADD);
	rneed=NeedParentheses(q,PRTY_ADD);
	r=(char *)MemNewN((int)StrLen(p)+lneed+1+(int)StrLen(q)+rneed+1);
	f=(char *) (lneed ? (rneed ? "(%s)*(%s)" : "(%s)*%s") :
		  (rneed ? "%s*(%s)" : "%s*%s"));
	sprintf(r,f,p,q);
      }
    }
  } else r=NULL;
  switch (freeptrs) {
    case 'l':
      MemFree(p);
      break;
    case 'r':
      MemFree(q);
      break;
    case 'b':
      MemFree(p);
      MemFree(q);
  }
  return r;
}

// f*g
Gradient operator*(const Gradient &f, const Gradient &g) {
  int differ;
  Gradient t;
  // construct expression f*g
  t.Expr[0]=Mult(f.Expr[f.Index],g.Expr[g.Index],0);
  differ=StrCmp(f.Expr[f.Index],g.Expr[g.Index]);
  switch (_Hod_) {
    case 3:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++)
	  for (k=j; k<_Niv_; k++) {
	    if (i==k)		// i==j==k
	      if (differ)
		desc="$1=$f3iii*$g0"    "$2=$f0*$g3iii"    "$1=$1+$2"
		     "$2=$f2ii*$g1i"    "$3=$f1i*$g2ii"    "$2=$2+$3"
					"$2=$c3*$2"
		     "$1=$1+$2";
	      else
		desc="$1=$c6*$f1i"      "$1=$1*$f2ii"
		     "$2=$c2*$f0"       "$2=$2*$f3iii"
		     "$1=$1+$2";
	    else if (i==j)	// i==j<k
	      if (differ)
		desc="$1=$f3iik*$g0"    "$2=$f0*$g3iik"    "$1=$1+$2"
		     "$2=$f2ii*$g1k"	"$1=$1+$2"
		     "$2=$f1k*$g2ii"     "$1=$1+$2"
		     "$2=$f2ik*$g1i"    "$3=$f1i*$g2ik"    "$2=$2+$3"
					"$2=$c2*$2"
		     "$1=$1+$2";
	      else
		desc="$1=$c4*$f1i"      "$1=$1*$f2ik"
		     "$2=$f1k*$f2ii"    "$3=$f0*$f3iik"    "$2=$2+$3"
					"$2=$c2*$2"
		     "$1=$1+$2";
	    else if (j==k)	// i<j==k
	      if (differ)
		desc="$1=$f3ijj*$g0"    "$2=$f0*$g3ijj"    "$1=$1+$2"
		     "$2=$f2jj*$g1i"	"$1=$1+$2"
		     "$2=$f1i*$g2jj"     "$1=$1+$2"
		     "$2=$f2ij*$g1j"    "$3=$f1j*$g2ij"    "$2=$2+$3"
					"$2=$c2*$2"
		     "$1=$1+$2";
	      else
		desc="$1=$c4*$f1j"      "$1=$1*$f2ij"
		     "$2=$f1i*$f2jj"    "$3=$f0*$f3ijj"    "$2=$2+$3"
					"$2=$c2*$2"
		     "$1=$1+$2";
	    else		// i<j<k
	      if (differ)
		desc="$1=$f3ijk*$g0"    "$2=$f0*$g3ijk"    "$1=$1+$2"
		     "$2=$f2ij*$g1k"    "$1=$1+$2"
		     "$2=$f2ik*$g1j"    "$1=$1+$2"
		     "$2=$f2jk*$g1i"    "$1=$1+$2"
		     "$2=$f1i*$g2jk"    "$1=$1+$2"
		     "$2=$f1j*$g2ik"    "$1=$1+$2"
		     "$2=$f1k*$g2ij"    "$1=$1+$2";
	      else
		desc="$1=$f1i*$f2jk"    "$2=$f1j*$f2ik"    "$1=$1+$2"
		     "$2=$f1k*$f2ij"    "$1=$1+$2"
		     "$2=$f0*$f3ijk"    "$1=$1+$2"
		     "$1=$c2*$1";
	    t.Der3[0][ind3(i,j,k)]=Interpret(f,g);
	  }
    case 2:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++) {
	  if (i==j)
	    if (differ)
	      desc="$1=$f2ii*$g0"     "$2=$f0*$g2ii"     "$1=$1+$2"
		   "$2=$c2*$f1i"      "$2=$2*$g1i"       "$1=$1+$2";
	    else
	      desc="$1=$f1i*$f1i"     "$2=$f0*$f2ii"     "$1=$1+$2"
		   "$1=$c2*$1"; 
	  else
	    if (differ)
	      desc="$1=$f2ij*$g0"     "$2=$f0*$g2ij"     "$1=$1+$2"
		   "$2=$f1i*$g1j"     "$1=$1+$2"
		   "$2=$f1j*$g1i"     "$1=$1+$2";
	    else
	      desc="$1=$f1i*$f1j"     "$2=$f0*$f2ij"     "$1=$1+$2"
		   "$1=$c2*$1";
	  t.Der2[0][ind2(i,j)]=Interpret(f,g);
	}
    case 1:
      // construct (f*g)'=f'*g+f*g'
      for (i=0; i<_Niv_; i++) {
	if (differ) desc="$1=$f1i*$g0"  "$2=$f0*$g1i"  "$1=$1+$2";
	else desc="$1=$c2*$f0"   "$1=$1*$f1i";
	t.Der1[0][i]=Interpret(f,g);
      }
  }
  return t;
}

// f*const
Gradient operator*(const Gradient &f, double c) {
  Gradient t(c);
  return f*t;
}

// const*f
Gradient operator*(double c, const Gradient &g) {
  Gradient t(c);
  return t*g;
}

// f*intconst
Gradient operator*(const Gradient &f, int ic) {
  Gradient t((double)ic);
  return f*t;
}

// intconst*f
Gradient operator*(int ic, const Gradient &g) {
  Gradient t((double)ic);
  return t*g;
}


//======= S u b t r a c t =======================

static char *Sub(char *p, char *q, char freeptrs) {
  char *r;
  char *f;
  int len,rneed;
  if (p) {
    len=(int)StrLen(p)+1;
    if (q) {
      rneed=NeedParentheses(q,PRTY_ADD);
      r=(char *)MemNewN(len+(int)StrLen(q)+rneed+1);
      f=(char *) (rneed ? "%s-(%s)" : "%s-%s");
      sprintf(r,f,p,q);
    } else {
      r=(char *)StringSave(p);
    }
  } else {
    if (q) {
      rneed=NeedParentheses(q,PRTY_ADD);
      r=(char *)MemNewN(1+(int)StrLen(q)+rneed+1);
      f=(char *) (rneed ? "-(%s)" : "-%s");
      sprintf(r,f,q);
    } else {
      r=NULL;
    }
  }
  switch (freeptrs) {
    case 'l':
      MemFree(p);
      break;
    case 'r':
      MemFree(q);
      break;
    case 'b':
      MemFree(p);
      MemFree(q);
  }
  return r;
}

// f-g
Gradient operator-(const Gradient &f, const Gradient &g) {
  Gradient t;
  // construct expression f-g
  t.Expr[0]=Sub(f.Expr[f.Index],g.Expr[g.Index],0);
  switch (_Hod_) {
    case 3:
      // construct (f-g)'''=f'''-g'''
      for (i=0; i<f.n3; i++)
	t.Der3[0][i]=Sub(f.Der3[f.Index][i],g.Der3[g.Index][i],0);
    case 2:
      // construct (f-g)''=f''-g''
      for (i=0; i<f.n2; i++)
	t.Der2[0][i]=Sub(f.Der2[f.Index][i],g.Der2[g.Index][i],0);
    case 1:
      // construct (f-g)'=f'-g'
      for (i=0; i<_Niv_; i++)
	t.Der1[0][i]=Sub(f.Der1[f.Index][i],g.Der1[g.Index][i],0);
  }
  return t;
}

// f-const
Gradient operator-(const Gradient &f, double c) {
  Gradient t(c);
  return f-t;
}

// -f
Gradient operator-(const Gradient &f) {
  return 0.0-f;
}

// const-f
Gradient operator-(double c, const Gradient &g) {
  Gradient t(c);
  return t-g;
}

// f-intconst
Gradient operator-(const Gradient &f, int ic) {
  Gradient t((double)ic);
  return f-t;
}

// intconst-f
Gradient operator-(int ic, const Gradient &g) {
  Gradient t((double)ic);
  return t-g;
}


//======= D i v i d e ===========================

static char *Div(char *p, char *q, char freeptrs) {
  char *r;
  char *f;
  int lneed,rneed;
  if (p && q) {
    if (EqualToOne(q)) r=(char *)StringSave(p);
    else {
      lneed=NeedParentheses(p,PRTY_ADD);
      rneed=NeedParentheses(q,PRTY_MULT);
      r=(char *)MemNewN((int)StrLen(p)+lneed+1+(int)StrLen(q)+rneed+1);
      f=(char *) (lneed ? (rneed ? "(%s)/(%s)" : "(%s)/%s") :
		(rneed ? "%s/(%s)" : "%s/%s"));
      sprintf(r,f,p,q);
    }
  } else r=NULL;
  switch (freeptrs) {
    case 'l':
      MemFree(p);
      break;
    case 'r':
      MemFree(q);
      break;
    case 'b':
      MemFree(p);
      MemFree(q);
  }
  return r;
}

// f/g
Gradient operator/(const Gradient &f, const Gradient &g) {
  Gradient t;
  // construct expression f/g
  t.Expr[0]=Div(f.Expr[f.Index],g.Expr[g.Index],0);
  switch (_Hod_) {
    case 3:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++)
	  for (k=j; k<_Niv_; k++) {
	    if (i==k)		// i==j==k
	      desc="$1=$g1i*$g1i"   "$1=$1*$g1i"   "$1=$1*$f0"   "$1=$1/$g0"
		   "$2=$f0*$g2ii"   "$3=$f1i*$g1i" "$2=$2+$3"    "$2=$2*$g1i"
				    "$1=$2-$1"     "$1=$c6*$1"   "$1=$1/$g0"
		   "$2=$f1i*$g2ii"  "$3=$f2ii*$g1i""$2=$2+$3"    "$2=$c3*$2"
				    "$3=$f0*$g3iii""$2=$2+$3"    "$1=$1-$2"
				    "$1=$1/$g0"
		   "$1=$f3iii+$1"   "$1=$1/$g0";
	    else if (i==j)	// i==j<k
	      desc="$1=$g1i*$g1i"   "$1=$1*$g1k"   "$1=$1*$f0"   "$1=$c6*$1"
				    "$1=$1/$g0"
		   "$2=$c4*$g1i"    "$2=$2*$g2ik"  "$3=$c2*$g1k" "$3=$3*$g2ii"
				    "$2=$2+$3"     "$2=$f0*$2"   "$1=$2-$1"
				    "$2=$c2*$f1k"  "$2=$2*$g1i"  "$3=$c4*$f1i"
				    "$3=$3*$g1k"   "$2=$2+$3"    "$2=$g1i*$2"
				    "$1=$1+$2"     "$1=$1/$g0"
		   "$2=$f0*$g3iik"  "$1=$1-$2"     "$2=$f2ii*$g1k""$1=$1-$2"
				    "$2=$f1k*$g2ii""$1=$1-$2"    "$2=$f1i*$g2ik"
				    "$3=$f2ik*$g1i""$2=$2+$3"    "$2=$c2*$2"
				    "$1-$1-$2"     "$1=$1/$g0"
		   "$1=$f3iik+$1"   "$1=$1/$g0";
	    else if (j==k)	// i<j==k
	      desc="$1=$g1j*$g1j"   "$1=$1*$g1i"   "$1=$1*$f0"   "$1=$c6*$1"
				    "$1=$1/$g0"
		   "$2=$c4*$g1j"    "$2=$2*$g2ij"  "$3=$c2*$g1i" "$3=$3*$g2jj"
				    "$2=$2+$3"     "$2=$f0*$2"   "$1=$2-$1"
				    "$2=$c2*$f1i"  "$2=$2*$g1j"  "$3=$c4*$f1j"
				    "$3=$3*$g1i"   "$2=$2+$3"    "$2=$g1j*$2"
				    "$1=$1+$2"     "$1=$1/$g0"
		   "$2=$f0*$g3ijj"  "$1=$1-$2"     "$2=$f2jj*$g1i""$1=$1-$2"
				    "$2=$f1i*$g2jj""$1=$1-$2"    "$2=$f1j*$g2ij"
				    "$3=$f2ij*$g1j""$2=$2+$3"    "$2=$c2*$2"
				    "$1-$1-$2"     "$1=$1/$g0"
		   "$1=$f3ijj+$1"   "$1=$1/$g0";
	    else		// i<j<k
	      desc="$1=$f0*$g1i"    "$1=$1*$g1j"   "$1=$1*$g1k"  "$1=$c6*$1"
				    "$1=$1/$g0"
		   "$2=$f0*$g2jk"   "$3=$f1k*$g1j" "$2=$2+$3"    "$4=$g1i*$2"
				    "$2=$f0*$g2ik" "$3=$f1i*$g1k""$2=$2+$3"
				    "$2=$g1j*$2"   "$4=$4+$2"    "$2=$f0*$g2ij"
				    "$3=$f1j*$g1i" "$2=$2+$3"    "$2=$g1k*$2"
				    "$4=$4+$2"     "$4=$c2*$4"   "$1=$4-$1"
				    "$1=$1/$g0"
		   "$2=$f0*$g3ijk"  "$3=$f1i*$g2jk""$2=$2+$3"    "$3=$f1j*$g2ik"
				    "$2=$2+$3"     "$3=$f1k*$g2ij""$2=$2+$3"
				    "$3=$g1i*$f2jk""$2=$2+$3"    "$3=$g1j*$f2ik"
				    "$2=$2+$3"     "$3=$g1k*$f2ij""$2=$2+$3"
				    "$1=$1-$2"     "$1=$1/$g0"   "$1=$f3ijk+$1"
				    "$1=$1/$g0";
	    t.Der3[0][ind3(i,j,k)]=Interpret(f,g);
	  }
    case 2:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++) {
	  if (i==j)
	    desc="$1=$g1i*$g1i"     "$1=$1*$f0"    "$1=$1*$c2"   "$1=$1/$g0"
		 "$2=$c2*$f1i"      "$2=$2*$g1i"   "$3=$f0*$g2ii""$2=$2+$3"
				    "$1=$1-$2"     "$1=$1/$g0"
		 "$1=$f2ii+$1"      "$1=$1/$g0";
	  else
	    desc="$1=$c2*$f0"       "$1=$1*$g1i"   "$1=$1*$g1j"  "$1=$1/$g0"
		 "$2=$f1i*$g1j"     "$3=$f1j*$g1i" "$2=$2+$3"    "$3=$f0*$g2ij"
				    "$2=$2+$3"     "$1=$1-$2"    "$1=$1/$g0"
		 "$1=$f2ij+$1"      "$1=$1/$g0";
	  t.Der2[0][ind2(i,j)]=Interpret(f,g);
	}
    case 1:
      // construct (f/g)'=(f'*g-f*g')/(g*g)=(f'-f*g'/g)/g
      for (i=0; i<_Niv_; i++) {
	desc="$1=$f0*$g1i"    "$1=$1/$g0"    "$1=$f1i-$1"  "$1=$1/$g0";
	t.Der1[0][i]=Interpret(f,g);
      }
  }
  return t;
}

// f/const
Gradient operator/(const Gradient &f, double c) {
  Gradient t(c);
  return f/t;
}

// const/f
Gradient operator/(double c, const Gradient &g) {
  Gradient t(c);
  return t/g;
}

// f/intconst
Gradient operator/(const Gradient &f, int ic) {
  return f/(double)ic;
}

// intconst/f
Gradient operator/(int ic, const Gradient &g) {
  return (double)ic/g;
}


//======= A s s i g n s =========================

// Important: for array elements in left-hand side
//	    of assignments (v[i]=f) we need
//	    access to v, not its copy v[i].
//	    Parent member provides a pointer to v.

// v=f
Gradient & Gradient::operator=(const Gradient &f) {
#if ED
  int i1,i2,i3,i4,i5;
  char name[100];
#endif	/* ED */
  char *p;
  char buf[50];
  // Output assignments and record where gradient is now
  switch (Kind) {
    case INDEP:
      fprintf(MStream_,"\n\n#error Assignment to independent variable is not allowed\n\n");
      break;
    case TEMP:
      // Output assignment operators
      p=f.Expr[f.Index];
#if ED
      switch (_Hod_) {
	case 4:
	case 5:
	  // Maple
	  fprintf(EDStream_,"%s:=%s:\n",Expr[Index],p ? p : "0");
	  return *this;
      }
#endif	/* ED */
      // Record that gradient is now in array associated with the temp
      fprintf(MStream_,"  %s=%s;\n",Expr[Index],p ? p : "0");
      switch (_Hod_) {
	case 3:
	  for (i=0; i<_Niv_; i++)
	    for (j=i; j<_Niv_; j++)
	      for (k=j; k<_Niv_; k++) {
		p=f.Der3[f.Index][ind3(i,j,k)];
		MemFree(Parent->Der3[Index][ind3(i,j,k)]);
		if (p) {
		  fprintf(MStream_,"  d3_%s[%i]=%s;\n",Expr[Index],ind3(i,j,k),p);
		  sprintf(buf,"d3_%s[%i]",Expr[Index],ind3(i,j,k));
		  Parent->Der3[Index][ind3(i,j,k)]=(char *)StringSave(buf);
		} else Parent->Der3[Index][ind3(i,j,k)]=NULL;
	      }
	case 2:
	  for (i=0; i<_Niv_; i++)
	    for (j=i; j<_Niv_; j++) {
	      p=f.Der2[f.Index][ind2(i,j)];
	      MemFree(Parent->Der2[Index][ind2(i,j)]);
	      if (p) {
		fprintf(MStream_,"  d2_%s[%i]=%s;\n",Expr[Index],ind2(i,j),p);
		sprintf(buf,"d2_%s[%i]",Expr[Index],ind2(i,j));
		Parent->Der2[Index][ind2(i,j)]=(char *)StringSave(buf);
	      } else Parent->Der2[Index][ind2(i,j)]=NULL;
	    }
	case 1:
	  for (i=0; i<_Niv_; i++) {
	    p=f.Der1[f.Index][i];
	    MemFree(Parent->Der1[Index][i]);
	    if (p) {
	      fprintf(MStream_,"  d1_%s[%i]=%s;\n",Expr[Index],i,p);
	      sprintf(buf,"d1_%s[%i]",Expr[Index],i);
	      Parent->Der1[Index][i]=(char *)StringSave(buf);
	    } else Parent->Der1[Index][i]=NULL;
	  }
      }
      fprintf(MStream_,"\n");
      break;
    case GRAD:
      // Output assignment operators
#if ED
      p=f.Expr[f.Index];
      switch (_Hod_) {
	case 4:
	  // Maple
	  fprintf(EDStream_,"\n");
	  // Specal processing if name is an array
	  i1=StrLen(Expr[Index])-1;
	  if (Expr[Index][i1]==']') {
	    sscanf(Expr[Index],"%[^[][%i]",name,&i2);
	    fprintf(EDStream_,"%s:=array(0..%i,0..%i):\n",
		name,i2,ind4(_Niv_,_Niv_,_Niv_,_Niv_));
	    sprintf(name,"%.*s,",i1,Expr[Index]);
	  } else {
	    sprintf(name,"%s[",Expr[Index]);
	    fprintf(EDStream_,"%s:=array(0..%i):\n",Expr[Index],ind4(_Niv_,_Niv_,_Niv_,_Niv_));
	  }
	  fprintf(EDStream_,"_expr_:=%s: _list_:=[]:\n\n",p ? p : "0");	// expr to differentiate
	  for (i1=0; i1<_Niv_; i1++) {
	    fprintf(EDStream_,"_d1_:=diff(_expr_,%s):\n",_IndepName[i1]);
	    for (i2=i1; i2<_Niv_; i2++) {
	      fprintf(EDStream_," _d2_:=diff(_d1_,%s):\n",_IndepName[i2]);
	      for (i3=i2; i3<_Niv_; i3++) {
		fprintf(EDStream_,"  _d3_:=diff(_d2_,%s):\n",_IndepName[i3]);
		for (i4=i3; i4<_Niv_; i4++) {
		  fprintf(EDStream_,"   _d4_:=diff(_d3_,%s):\n",_IndepName[i4]);
		  fprintf(EDStream_,"   if _d4_<>0.0 then _list_:=[op(_list_),%s%i]=_d4_]: fi:\n",
			name,ind4(i1,i2,i3,i4));
		}
	      }
	    }
	  }
	  fprintf(EDStream_,"if _list_<>[] then C(_list_,optimized): fi:\n");
	  return *this;
	case 5:
	  // Maple
	  fprintf(EDStream_,"\n");
	  // Specal processing if name is an array
	  i1=StrLen(Expr[Index])-1;
	  if (Expr[Index][i1]==']') {
	    sscanf(Expr[Index],"%[^[][%i]",name,&i2);
	    fprintf(EDStream_,"%s:=array(0..%i,0..%i):\n",
		name,i2,ind5(_Niv_,_Niv_,_Niv_,_Niv_,_Niv_));
	    sprintf(name,"%.*s,",i1,Expr[Index]);
	  } else {
	    sprintf(name,"%s[",Expr[Index]);
	    fprintf(EDStream_,"%s:=array(0..%i):\n",Expr[Index],ind5(_Niv_,_Niv_,_Niv_,_Niv_,_Niv_));
	  }
	  fprintf(EDStream_,"_expr_:=%s: _list_:=[]:\n\n",p ? p : "0");	// expr to differentiate
	  for (i1=0; i1<_Niv_; i1++) {
	    fprintf(EDStream_,"_d1_:=diff(_expr_,%s):\n",_IndepName[i1]);
	    for (i2=i1; i2<_Niv_; i2++) {
	      fprintf(EDStream_," _d2_:=diff(_d1_,%s):\n",_IndepName[i2]);
	      for (i3=i2; i3<_Niv_; i3++) {
		fprintf(EDStream_,"  _d3_:=diff(_d2_,%s):\n",_IndepName[i3]);
		for (i4=i3; i4<_Niv_; i4++) {
		  fprintf(EDStream_,"   _d4_:=diff(_d3_,%s):\n",_IndepName[i4]);
		  for (i5=i4; i5<_Niv_; i5++) {
		    fprintf(EDStream_,"    _d5_:=diff(_d4_,%s):\n",_IndepName[i5]);
		    fprintf(EDStream_,"    if _d5_<>0.0 then _list_:=[op(_list_),%s%i]=_d5_]: fi:\n",
			name,ind5(i1,i2,i3,i4,i5));
		  }
		}
	      }
	    }
	  }
	  fprintf(EDStream_,"if _list_<>[] then C(_list_,optimized): fi:\n");
	  return *this;
      }
#endif	/* ED */
      // Gradient: just copy
      // NB: Don't assign zero because all the gradient vars are initially set to 0
      switch (_Hod_) {
	case 3:
	  for (i=0; i<_Niv_; i++)
	    for (j=i; j<_Niv_; j++)
	      for (k=j; k<_Niv_; k++) {
		p=f.Der3[f.Index][ind3(i,j,k)];
		if (p)
		  fprintf(MStream_,"  %s[%i]=%s; /* /d%id%id%i */\n",Expr[Index],ind3(i,j,k),p,i,j,k);
		MemFree(Parent->Der3[Index][ind3(i,j,k)]);
		Parent->Der3[Index][ind3(i,j,k)]=(char *)StringSave(p);
	      }
	  break;
	case 2:
	  for (i=0; i<_Niv_; i++)
	    for (j=i; j<_Niv_; j++) {
	      p=f.Der2[f.Index][ind2(i,j)];
	      if (p)
		fprintf(MStream_,"  %s[%i]=%s; /* /d%id%i */\n",Expr[Index],ind2(i,j),p,i,j);
	      MemFree(Parent->Der2[Index][ind2(i,j)]);
	      Parent->Der2[Index][ind2(i,j)]=(char *)StringSave(p);
	    }
	  break;
	case 1:
	  for (i=0; i<_Niv_; i++) {
	    p=f.Der1[f.Index][i];
	    if (p)
	      fprintf(MStream_,"  %s[%i]=%s;\n",Expr[Index],i,p);
	    MemFree(Parent->Der1[Index][i]);
	    Parent->Der1[Index][i]=(char *)StringSave(p);
	  }
	  fprintf(MStream_,"\n");
      }
      break;
    case DEF:
      // don't output assign statements
      switch (_Hod_) {
	case 3:
	  for (i=0; i<_Niv_; i++)
	    for (j=i; j<_Niv_; j++)
	      for (k=j; k<_Niv_; k++) {
		MemFree(Der3[Index][ind3(i,j,k)]);
		Der3[Index][ind3(i,j,k)]=(char *)StringSave(f.Der3[f.Index][ind3(i,j,k)]);
	      }
	  break;
	case 2:
	  for (i=0; i<_Niv_; i++)
	    for (j=i; j<_Niv_; j++) {
	      MemFree(Der2[Index][ind2(i,j)]);
	      Der2[Index][ind2(i,j)]=(char *)StringSave(f.Der2[f.Index][ind2(i,j)]);
	    }
	  break;
	case 1:
	  for (i=0; i<_Niv_; i++) {
	    MemFree(Der1[Index][i]);
	    Der1[Index][i]=(char *)StringSave(f.Der1[f.Index][i]);
	  }
      }
      MemFree(Expr[Index]);
      Expr[Index]=(char *)StringSave(f.Expr[f.Index]);
      break;
    default:;
  }
  ResetIndex();
  return *this;
}

// v=const
Gradient & Gradient::operator=(double c) {
  Gradient t(c);
  *this=t;
  return *this;
}

// v=intconst
Gradient & Gradient::operator=(int ic) {
  *this=(double)ic;
  return *this;
}

// v+=f;
Gradient & Gradient::operator+=(const Gradient &f) {
  *this=*this+f;
  return *this;
}

// v+=const
Gradient & Gradient::operator+=(double c) {
  *this=*this+c;
  return *this;
}

// v+=intconst
Gradient & Gradient::operator+=(int ic) {
  *this=*this+ic;
  return *this;
}

// v-=f
Gradient & Gradient::operator-=(const Gradient &f) {
  *this=*this-f;
  return *this;
}

// v-=const
Gradient & Gradient::operator-=(double c) {
  *this=*this-c;
  return *this;
}

// v-=intconst
Gradient & Gradient::operator-=(int ic) {
  *this=*this-ic;
  return *this;
}

// v*=f
Gradient & Gradient::operator*=(const Gradient &f) {
  *this=*this*f;
  return *this;
}

// v*=const
Gradient & Gradient::operator*=(double c) {
  *this=*this*c;
  return *this;
}

// v*=intconst
Gradient & Gradient::operator*=(int ic) {
  *this=*this*ic;
  return *this;
}


// v/=f
Gradient & Gradient::operator/=(const Gradient &f) {
  *this=*this/f;
  return *this;
}

// v/=const
Gradient & Gradient::operator/=(double c) {
  *this=*this/c;
  return *this;
}

// v/=intconst
Gradient & Gradient::operator/=(int ic) {
  *this=*this/ic;
  return *this;
}


//======= F u n c t i o n s =====================

// Creates texts of all partial derivatives of stdfun(arg) assuming that
//  d1=dstdfun/darg, d2=d2stdfun/darg2, d3=d3stdfun/darg3
static void DerText(Gradient &res, const Gradient &arg, char *d1, char *d2, char *d3) {
  char *p;
  char *q;
  char *w;
  switch (_Hod_) {
    case 3:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++)
	  for (k=j; k<_Niv_; k++) {
	    desc="$1=$f1i*$f1j"   "$1=$1*$f1k";	//arg'x*arg'y*arg'z
	    p=Interpret(arg,arg);
	    p=Mult(p,d3,'l');
	    q=Mult(GetDer3T_(arg,i,j,k),d1,0);
	    p=Add(p,q,'b');			//D'''*arg'x*arg'y*arg'z+D'*arg'''xyz
	    if (i==k) {		// i==j==k
	      desc="$1=$c3*$f1i"   "$1=$1*$f2ii";
	      q=Interpret(arg,arg);
	    } else if (i==j) {	// i==j<k
	      desc="$1=$c2*$f1i"   "$1=$1*$f2ik"  "$2=$f1k*$f2ii" "$1=$1+$2";
	      q=Interpret(arg,arg);
	    } else if (j==k) {	// i<j==k
	      desc="$1=$c2*$f1j"   "$1=$1*$f2ij"  "$2=$f1i*$f2jj" "$1=$1+$2";
	      q=Interpret(arg,arg);
	    } else {		// i<j<k
	      desc="$1=$f1i*$f2jk"  "$2=$f1j*$f2ik" "$1=$1+$2"
				    "$2=$f1k*$f2ij" "$1=$1+$2";
	      q=Interpret(arg,arg);
	    }
	    q=Mult(q,d2,'l');
	    res.SetDer3(i,j,k,w=Add(p,q,'b'));
	    MemFree(w);
	  }
    case 2:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++) {
	  desc="$1=$f1i*$f1j";
	  p=Interpret(arg,arg);
	  p=Mult(p,d2,'l');
	  q=Mult(GetDer2T_(arg,i,j),d1,0);
	  res.SetDer2(i,j,w=Add(p,q,'b'));
	  MemFree(w);
	}
    case 1:
      for (i=0; i<_Niv_; i++) {
	res.SetDer1(i,w=Mult(d1,GetDer1T_(arg,i),0));
	MemFree(w);
      }
  }
}

static char val_buf[20],grad_buf[20];

static void Arrays(void) {
  int indx;
  indx=sfGetIndex();
  sprintf(val_buf,"sv__[%i]",indx);
  sprintf(grad_buf,"sd__[%i]",indx);
}

static void MultAll(Gradient &t, char *mult) {
  char *p;
  switch (_Hod_) {
    case 3:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++)
	  for (k=j; k<_Niv_; k++) {
	    t.SetDer3(i,j,k,p=Mult(mult,GetDer3T_(t,i,j,k),0));	// 0, not 'r'!
	    MemFree(p);
	  }
    case 2:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++) {
	  t.SetDer2(i,j,p=Mult(mult,GetDer2T_(t,i,j),0));	// 0, not 'r'!
	  MemFree(p);
	}
    case 1:
      for (i=0; i<_Niv_; i++) {
	t.SetDer1(i,p=Mult(mult,GetDer1T_(t,i),0));		// 0, not 'r'!
	MemFree(p);
      }
  }
}

Gradient fabs(const Gradient &f) {
  Gradient t;
  char *d1;
  int l=(int)StringLen(f.Expr[f.Index])+1;
  t.Expr[0]=(char *)MemNewN(l+6);
  sprintf(t.Expr[0],"fabs(%s)",f.Expr[f.Index]);
#if ED
  if (_Hod_>3) return t;
#endif	/* ED */
  d1=(char *)MemNewN(l+14);
  sprintf(d1,"(%s>=0 ? 1 : -1)",f.Expr[f.Index]);
  DerText(t,f,d1,NULL,NULL);
  MemFree(d1);
  return t;
}

Gradient exp(const Gradient &f) {
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+5);
    sprintf(t.Expr[0],"exp(%s)",f.Expr[f.Index]);
    return t;
  }
#endif	/* ED */
  Arrays();		// assign "sv__[]" to val_buf
  fprintf(MStream_,"  %s=exp(%s);\n",val_buf,f.Expr[f.Index]);
  t.Expr[0]=(char *)StrSave(val_buf);
  DerText(t,f,"1","1","1");
  MultAll(t,val_buf);
  return t;
}

Gradient log(const Gradient &f) {
  Gradient t;
  char *d1;
  char *d2;
  char *d3;
  t.Expr[0]=(char *)MemNewN((int)StringLen(f.Expr[f.Index])+5+1);
  sprintf(t.Expr[0],"log(%s)",f.Expr[f.Index]);
#if ED
  if (_Hod_>3) return t;
#endif	/* ED */
  desc="$1=$c1/$f0";			// log'(u)=1/u
  d1=Interpret(f,f);
  desc="$1=$f0*$f0" "$1=$c1/$1" "$1=$c0-$1";	// log''(u)=-1/(u*u)
  d2=Interpret(f,f);
  desc="$1=$f0*$f0" "$1=$1*$f0" "$1=$c2/$1";	// log'''(u)=2/(u*u*u)
  d3=Interpret(f,f);
  DerText(t,f,d1,d2,d3);
  MemFree(d1);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient sqrt(const Gradient &f) {
  char *d1;
  char *d2;
  char *d3;
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+6);
    sprintf(t.Expr[0],"sqrt(%s)",f.Expr[f.Index]);
    return t;
  }
#endif	/* ED */
  Arrays();		// assign "sv__[]" to val_buf and "sd__[]" to grad_buf
  fprintf(MStream_,"  %s=sqrt(%s);\n",val_buf,f.Expr[f.Index]);
  t.Expr[0]=(char *)StrSave(val_buf);
  fprintf(MStream_,"  %s=1/(2*%s);\n",grad_buf,val_buf);	//mult is 1/(2sqrt(u))
  d1="1";							// sqrt'(u)=1*mult
  desc="$1=$c2*$f0"  "$1=$c1/$1"   "$1=$c0-$1";	// sqrt''(u)=-1/(2u)*mult
  d2=Interpret(f,t);
  desc="$1=$f0*$f0"  "$1=$1*$c4"   "$1=$c3/$1";	// sqrt'''(u)=3/(4*u*u)*mult
  d3=Interpret(f,f);
  DerText(t,f,d1,d2,d3);
  MultAll(t,grad_buf);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient pow(const Gradient &f, double p) {
  char *d1;
  char *d2;
  char *d3;
  char *mult;
  char buf[20];
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+60);
    sprintf(t.Expr[0],"pow(%s,%g)",f.Expr[f.Index],p);
    return t;
  }
#endif	/* ED */
  if (p==2) {
    desc="$1=$f0*$f0";
    t.Expr[0]=Interpret(f,f);
  } else {
    sprintf(buf,"%g",p);
    t.Expr[0]=(char *)MemNewN((int)StringLen(f.Expr[f.Index])+(int)StringLen(buf)+6+1);
    sprintf(t.Expr[0],"pow(%s,%s)",f.Expr[f.Index],buf);
  }
  d1="1";
  if (p==1) {
    mult=d2=d3=NULL;
  } else
    if (p==2) {
      mult=Mult("2",f.Expr[f.Index],0);
      d2=Div("1",f.Expr[f.Index],0);
      d3=NULL;
    } else {
      Arrays();		// send "sd__[]" to grad_buf
      fprintf(MStream_,"  %s=%g*pow(%s,%g);\n",grad_buf,p,f.Expr[f.Index],p-1);
      mult=(char *)StrSave(grad_buf);
      sprintf(buf,"%g",p-1);
      d2=Div(buf,f.Expr[f.Index],0);
      sprintf(buf,"%g",(p-1)*(p-2));
      d3=Div(buf,Mult(f.Expr[f.Index],f.Expr[f.Index],0),'r');
    }
  DerText(t,f,d1,d2,d3);
  if (p!=1) MultAll(t,mult);
  MemFree(mult);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient pow(const Gradient &f, int ic) {
  return pow(f,(double)ic);
}

Gradient pow(const Gradient &f, const Gradient &g) {
  return exp(g*log(f));
}

Gradient sin(const Gradient &f) {
  char *d1;
  char *d2;
  char *d3;
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+5);
    sprintf(t.Expr[0],"sin(%s)",f.Expr[f.Index]);
    return t;
  }
#endif	/* ED */
  Arrays();	// "sv__[]" to val_buf and "sd__[]" to grad_buf
  fprintf(MStream_,"  %s=sin(%s);\n",val_buf,f.Expr[f.Index]);
  t.Expr[0]=(char *)StrSave(val_buf);
  fprintf(MStream_,"  %s=cos(%s);\n",grad_buf,f.Expr[f.Index]);
  d1=grad_buf;			// sin'=cos
  d2=Sub(NULL,val_buf,0);	// sin''=-sin
  d3=Sub(NULL,grad_buf,0);	// sin'''=-cos
  DerText(t,f,d1,d2,d3);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient cos(const Gradient &f) {
  char *d1;
  char *d2;
  char *d3;
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+5);
    sprintf(t.Expr[0],"cos(%s)",f.Expr[f.Index]);
    return t;
  }
#endif	/* ED */
  Arrays();	// "sv__[]" to val_buf and "sd__[]" to grad_buf
  fprintf(MStream_,"  %s=cos(%s);\n",val_buf,f.Expr[f.Index]);
  t.Expr[0]=(char *)StrSave(val_buf);
  fprintf(MStream_,"  %s=-sin(%s);\n",grad_buf,f.Expr[f.Index]);
  d1=grad_buf;			// cos'=-sin
  d2=Sub(NULL,val_buf,0);	// cos''=-cos
  d3=Sub(NULL,grad_buf,0);	// cos'''=sin
  DerText(t,f,d1,d2,d3);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient tan(const Gradient &f) {
  char *d1;
  char *d2;
  char *d3;
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+5);
    sprintf(t.Expr[0],"tan(%s)",f.Expr[f.Index]);
    return t;
  }
#endif	/* ED */
  Arrays();	// "sv__[]" to val_buf and "sd__[]" to grad_buf
  fprintf(MStream_,"  %s=tan(%s);\n",val_buf,f.Expr[f.Index]);
  t.Expr[0]=(char *)StrSave(val_buf);
  fprintf(MStream_,"  %s=1+tan(%s)*tan(%s);\n",grad_buf,f.Expr[f.Index],f.Expr[f.Index]);
					// mult is 1+tan*tan
  d1="1";				// tan'=1*mult
  d2=Mult("2",val_buf,0);		// tan''=(2*tan)*mult
  d3=Sub(Mult("6",grad_buf,0),"4",'l');	// tan'''=(6*mult-4)*mult
  DerText(t,f,d1,d2,d3);
  MultAll(t,grad_buf);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient atan(const Gradient &f) {
  char *d1;
  char *d2;
  char *d3;
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+8);
    sprintf(t.Expr[0],"arctan(%s)",f.Expr[f.Index]);
    return t;
  }
#endif	/* ED */
  Arrays();	// "sd__[]" to grad_buf
  t.Expr[0]=(char *)MemNewN((int)StringLen(f.Expr[f.Index])+6+1);
  sprintf(t.Expr[0],"atan(%s)",f.Expr[f.Index]);
  fprintf(MStream_,"  %s=1/(1+(%s)*(%s));\n",grad_buf,f.Expr[f.Index],f.Expr[f.Index]);
					// mult is 1/(1+u*u)
  d1="1";				// atan'(u)=1*mult
  d2=Mult("-2",Mult(f.Expr[f.Index],grad_buf,0),'r');	// atan''(u)=(-2*u*mult)*mult
  d3=Mult(grad_buf,Sub("6",Mult("8",grad_buf,0),'r'),'r');	// atan'''(u)=(mult(6-8mult))mult
  DerText(t,f,d1,d2,d3);
  MultAll(t,grad_buf);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient asin(const Gradient &f) {
  char *d1;
  char *d2;
  char *d3;
  Gradient t;
#if ED
  if (_Hod_>3) {
    int l=(int)StringLen(f.Expr[f.Index])+1;
    t.Expr[0]=(char *)MemNewN(l+8);
    sprintf(t.Expr[0],"arcsin(%s)",f.Expr[f.Index]);
    return t;
  }
#endif	/* ED */
  Arrays();	// "sv__[]" to val_buf and "sd__[]" to grad_buf
  t.Expr[0]=(char *)MemNewN((int)StringLen(f.Expr[f.Index])+6+1);
  sprintf(t.Expr[0],"asin(%s)",f.Expr[f.Index]);
  fprintf(MStream_,"  %s=1-(%s)*(%s);\n",val_buf,f.Expr[f.Index],f.Expr[f.Index]);
  fprintf(MStream_,"  %s=1/sqrt(%s);\n",grad_buf,val_buf);
					// mult is 1/sqrt(1-u*u)
  d1="1";				// asin'(u)=1*mult
  d2=Div(f.Expr[f.Index],val_buf,0);	// atan''(u)=u/(1-u*u)*mult
  d3=Sub("3",Mult("2",val_buf,0),'r');
  d3=Div(Div(d3,val_buf,'l'),val_buf,'l');	// asin'''(u)=(3-2(1-u*u)/(1-u*u)*(1-u*u))*mult
  DerText(t,f,d1,d2,d3);
  MultAll(t,grad_buf);
  MemFree(d2);
  MemFree(d3);
  return t;
}

Gradient acos(const Gradient &f) {
  return -asin(f);
}

// iff(c,f,g) is equivalent to  c ? f : g
Gradient iff(const Gradient &c, const Gradient &f, const Gradient &g) {
  Gradient t;
  char *fp;
  char *gp;
  char *p;
  int len;
  len=(int)StringLen(c.Expr[c.Index])+1;
  t.Expr[0]=(char *)MemNewN((int)(len+StringLen(f.Expr[f.Index])+StringLen(g.Expr[g.Index])+7));
  sprintf(t.Expr[0],"iff(%s,%s,%s)",c.Expr[c.Index],f.Expr[f.Index],g.Expr[g.Index]);
  switch (_Hod_) {
    case 3:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++)
	  for (k=j; k<_Niv_; k++) {
	    fp=GetDer3T_(f,i,j,k);
	    gp=GetDer3T_(g,i,j,k);
	    if (fp || gp) {
	      p=(char *)MemNewN((int)(len+StringLen(fp)+StringLen(gp)+7));
	      sprintf(p,"iff(%s,%s,%s)",c.Expr[c.Index],fp ? fp : "0",gp ? gp : "0");
	      t.SetDer3(i,j,k,p);
	    }
	  }
    case 2:
      for (i=0; i<_Niv_; i++)
	for (j=i; j<_Niv_; j++) {
	  fp=GetDer2T_(f,i,j);
	  gp=GetDer2T_(g,i,j);
	  if (fp || gp) {
	    p=(char *)MemNewN((int)(len+StringLen(fp)+StringLen(gp)+7));
	    sprintf(p,"iff(%s,%s,%s)",c.Expr[c.Index],fp ? fp : "0",gp ? gp : "0");
	    t.SetDer2(i,j,p);
	  }
	}
    case 1:
      for (i=0; i<_Niv_; i++) {
	fp=GetDer1T_(f,i);
	gp=GetDer1T_(g,i);
	if (fp || gp) {
	  p=(char *)MemNewN((int)(len+StringLen(fp)+StringLen(gp)+7));
	  sprintf(p,"iff(%s,%s,%s)",c.Expr[c.Index],GetDer1_(f,i),GetDer1_(g,i));
	  t.SetDer1(i,p);
	}
      }
  }
  return t;
}

//======= R E L A T I O N S    ( I F F ) ========

static char *rel(char *p, char *q, char *op, int prty) {
  char *r;
  char *f;
  int len,lneed,rneed;
  len=(int)(StringLen(p)+StrLen(op)+StringLen(q)+1);
  lneed=NeedParentheses(p,prty);
  rneed=NeedParentheses(q,prty);
  len+=lneed+rneed;
  r=(char *)MemNewN(len+lneed+rneed);
  f=(char *) (lneed ? (rneed ? "(%s)%s(%s)" : "(%s)%s%s") :
	    (rneed ? "%s%s(%s)" : "%s%s%s"));
  sprintf(r,f,p,op,q);
  return r;
}

// p <rel> q   <rel> is < <= > >= 
#define Rel(p,q,op) rel(p,q,op,PRTY_REL)


// p <equ> q   <equ> is == !=
#define Equ(p,q,op) rel(p,q,op,PRTY_EQU)

// f>g
Gradient operator>(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=Rel(f.Expr[f.Index],g.Expr[g.Index],">");
  return t;
}

// f>const
Gradient operator>(const Gradient &f, double c) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(f.Expr[f.Index],b,">");
  return t;
}

// f>intconst
Gradient operator>(const Gradient &f, int ic) {
  return f>(double)ic;
}

// const>f
Gradient operator>(double c, const Gradient &f) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(b,f.Expr[f.Index],">");
  return t;
}

// intconst>f
Gradient operator>(int ic, const Gradient &f) {
  return (double)ic>f;
}

// f>=g
Gradient operator>=(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=Rel(f.Expr[f.Index],g.Expr[g.Index],">=");
  return t;
}

// f>=const
Gradient operator>=(const Gradient &f, double c) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(f.Expr[f.Index],b,">=");
  return t;
}

// f>=intconst
Gradient operator>=(const Gradient &f, int ic) {
  return f>=(double)ic;
}

// const>=f
Gradient operator>=(double c, const Gradient &f) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(b,f.Expr[f.Index],">=");
  return t;
}

// intconst>=f
Gradient operator>=(int ic, const Gradient &f) {
  return (double)ic>=f;
}

// f<g
Gradient operator<(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=Rel(f.Expr[f.Index],g.Expr[g.Index],"<");
  return t;
}

// f<const 
Gradient operator<(const Gradient &f, double c) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(f.Expr[f.Index],b,"<");
  return t;
}

// f<intconst
Gradient operator<(const Gradient &f, int ic) {
  return f<(double)ic;
}

// const<f
Gradient operator<(double c, const Gradient &f) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(b,f.Expr[f.Index],"<");
  return t;
}

// intconst<f
Gradient operator<(int ic, const Gradient &f) {
  return (double)ic<f;
}


// f<=g
Gradient operator<=(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=Rel(f.Expr[f.Index],g.Expr[g.Index],"<=");
  return t;
}

// f<=const
Gradient operator<=(const Gradient &f, double c) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(f.Expr[f.Index],b,"<=");
  return t;
}

// f<=intconst
Gradient operator<=(const Gradient &f, int ic) {
  return f<=(double)ic;
}

// const <=f
Gradient operator<=(double c, const Gradient &f) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Rel(b,f.Expr[f.Index],"<=");
  return t;
}

// intconst<=f
Gradient operator<=(int ic, const Gradient &f) {
  return (double)ic<=f;
}


// f==g
Gradient operator==(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=Equ(f.Expr[f.Index],g.Expr[g.Index],"==");
  return t;
}

// f==const
Gradient operator==(const Gradient &f, double c) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Equ(f.Expr[f.Index],b,"==");
  return t;
}

// f==intconst
Gradient operator==(const Gradient &f, int ic) {
  return f==(double)ic;
}

// const==f
Gradient operator==(double c, const Gradient &f) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Equ(b,f.Expr[f.Index],"==");
  return t;
}

// intconst==f
Gradient operator==(int ic, const Gradient &f) {
  return (double)ic==f;
}


// f!=g
Gradient operator!=(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=Equ(f.Expr[f.Index],g.Expr[g.Index],"!=");
  return t;
}

// f!=const
Gradient operator!=(const Gradient &f, double c) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Equ(f.Expr[f.Index],b,"!=");
  return t;
}

// f!=intconst
Gradient operator!=(const Gradient &f, int ic) {
  return f!=(double)ic;
}

// const!=f
Gradient operator!=(double c, const Gradient &f) {
  Gradient t;
  char b[20];
  sprintf(b,"%g",c);
  t.Expr[0]=Equ(b,f.Expr[f.Index],"!=");
  return t;
}

// intconst!=f
Gradient operator!=(int ic, const Gradient &f) {
  return (double)ic!=f;
}

//======= B O O L E A N S =======================

// !f
Gradient operator!(const Gradient &f) {
  Gradient t;
  int need;
  need=NeedParentheses(f.Expr[f.Index],PRTY_UNARY);
  t.Expr[0]=(char *)MemNewN((int)(1+StrLen(f.Expr[f.Index])+need+1));
  sprintf(t.Expr[0],need ? "!(%s)" : "!%s",f.Expr[f.Index]);
  return t;
}

// p&&q
#define BoolAnd(p,q) rel(p,q,"&&",PRTY_AND)

// f&&g
Gradient operator&&(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=BoolAnd(f.Expr[f.Index],g.Expr[g.Index]);
  return t;
}

// p||q
#define BoolOr(p,q) rel(p,q,"||",PRTY_OR)

// f||g
Gradient operator||(const Gradient &f, const Gradient &g) {
  Gradient t;
  t.Expr[0]=BoolOr(f.Expr[f.Index],g.Expr[g.Index]);
  return t;
}

//======= A c c e s s   t o   c o m p o n e n t s

// caller must not change or freeptrs memory
char *GetExpr_(const Gradient &f) {
  char *p;
  p=f.Expr[f.Index];
  return (char *) (p ? p : "0");
}

char *GetExprT_(const Gradient &f) {
  return f.Expr[f.Index];
}

void Gradient::SetExpr(char *newtext) {
  MemFree(this->Expr[this->Index]);
  this->Expr[this->Index]=(char *)StringSave(newtext);
}

// caller must not change or freeptrs memory
char *GetDer1_(const Gradient &f, int ic) {
  char *p;
  p=f.Der1[f.Index][ic];
  return (char *) (p ? p : "0");
}

char *GetDer1T_(const Gradient &f, int ic) {
  return f.Der1[f.Index][ic];
}

void Gradient::SetDer1(int ic, char *newtext) {
  MemFree(this->Der1[this->Index][ic]);
  this->Der1[this->Index][ic]=(char *)StringSave(newtext);
}

// caller must not change or freeptrs memory
char *GetDer2_(const Gradient &f, int ic, int jc) {
  char *p;
  p=f.Der2[f.Index][ind2(ic,jc)];
  return (char *) (p ? p : "0");
}

char *GetDer2T_(const Gradient &f, int ic, int jc) {
  return f.Der2[f.Index][ind2(ic,jc)];
}

void Gradient::SetDer2(int ic, int jc, char *newtext) {
  int indx=ind2(ic,jc);
  MemFree(this->Der2[this->Index][indx]);
  this->Der2[this->Index][indx]=(char *)StringSave(newtext);
}

// caller must not change or freeptrs memory
char *GetDer3_(const Gradient &f, int ic, int jc, int kc) {
  char *p;
  p=f.Der3[f.Index][ind3(ic,jc,kc)];
  return (char *) (p ? p : "0");
}

char *GetDer3T_(const Gradient &f, int ic, int jc, int kc) {
  return f.Der3[f.Index][ind3(ic,jc,kc)];
}

void Gradient::SetDer3(int ic, int jc, int kc, char *newtext) {
  int indx=ind3(ic,jc,kc);
  MemFree(this->Der3[this->Index][indx]);
  this->Der3[this->Index][indx]=(char *)StringSave(newtext);
}


// Interpret a sequence of assignments
// Left-hand side: $1, $2, $3, $4
// Right-hand side: binary operation (+,-,*,/), operands may be
//		  $1, $2, $3, $4, $cn (n is any integer, n>0), $f0, $g0,
//		  $f1d, $g1d, $f2dd, $g2dd, $f3ddd, $g3ddd
//		  where d may by any letter i,j,k

static char Letter[]={0,'r','l','b'};

static char *Interpret(const Gradient &f, const Gradient &g) {
  char *tmp[4];
  char *op[2];
  char *(*Oper)(char*,char*,char);
  int o,ii,jj,pos;
  char freeptrs[2];
  char buf[10];
  while (*desc) {
    ii=desc[1]-'1';
    desc+=4;
    for (jj=0; 1; jj++) {
      switch (*desc) {
	case 'f':
	case 'g':
	  o=desc[1]-'0';
	  switch (o) {
	    case 0:
	      op[jj]=GetExprT_(*desc=='f' ? f : g);
	      break;
	    case 1:
	      op[jj]=GetDer1T_(*desc=='f' ? f : g,IndeX[desc[2]-'i']);
	      break;
	    case 2:
	      op[jj]=GetDer2T_(*desc=='f' ? f : g,IndeX[desc[2]-'i'],IndeX[desc[3]-'i']);
	      break;
	    case 3:
	      op[jj]=GetDer3T_(*desc=='f' ? f : g,IndeX[desc[2]-'i'],IndeX[desc[3]-'i'],IndeX[desc[4]-'i']);
	      break;
	  }
	  desc+=o+2;
	  freeptrs[jj]=0;
	  break;
	case 'c':
	  sscanf(desc+1,"%[0-9]%n",buf,&pos);
	  if (StrCmp(buf,"0")) op[jj]=(char *)StringSave(buf);
	  else op[jj]=0;
	  desc+=pos+1;
	  freeptrs[jj]=1;
	  break;
	default:	/* must be '1', '2', '3', or '4' */
	  op[jj]=tmp[*desc-'1'];
	  desc+=1;
	  freeptrs[jj]=1;
      }
      if (jj) break;
      else {
	switch (*desc) {
	  case '+': Oper=Add; break;
	  case '-': Oper=Sub; break;
	  case '*': Oper=Mult; break;
	  case '/': Oper=Div; break;
	}
	desc+=2;
      }
    }
    tmp[ii]=Oper(op[0],op[1],Letter[2*freeptrs[0]+freeptrs[1]]);
  }
  return tmp[ii];
}


//***********************************//
// C++ run-time processing functions //
//***********************************//


static FILE *SaveStream[MAXDEP_];

// Initializes run-time support
void rtInit_(FILE *m, FILE *f) {
  MStream_=m;
  FStream_=f;
  FuncLevel_=0;
#if ED
  if (_Hod_<=3) /* skip for External Differentiator (i.e. Maple) */
#endif
  sf_set[0].Max=uf_set[0].Max=1;
  Memset(sf_set[0].CurSet,0,sizeof(Set));
  Memset(sf_set[0].OrSet,0,sizeof(Set));
  Memset(uf_set[0].CurSet,0,sizeof(Set));
  Memset(uf_set[0].OrSet,0,sizeof(Set));
}

// Terminates run-time support
void rtTerm_(void) {
#if ED
  {
    int ii;
    if (_IndepNum) {
      for (ii=0; ii<_IndepNum; ii++) {
	MemFree(_IndepName[ii]);
      }
      MemFree(_IndepName);
      _IndepName=NULL;
      _IndepNum=0;
    }
  }
#endif	/* ED */
}

#if ED
// Called before (PreProc==1) or after (PreProc==0) generation of
// the text for External Differentiator
void rtExtDiff_(char *mio, int PreProc) {
  switch (PreProc) {
    case 1:	 // before output
      // Open
      EDStream_=fopen(mio,"wt");
      if (!EDStream_) {
	fprintf(MStream_,"#error Cannot open file '%s'\n",mio);
	return;
      }
      // Maple-specific
      fprintf(EDStream_,"readlib(C):\n");
      fprintf(EDStream_,"pow:=(x,d)->x^d:\n");
      fprintf(EDStream_,"fabs:=(x)->abs(x):\n");
      fprintf(EDStream_,"iff:=proc(c,e1,e2) if c<>0 then e1 else e2  fi end:\n");
      break;
    case 0:	// after output
      if (!EDStream_) return;
      // Maple-specific
      fprintf(EDStream_,"quit:\n");
      // Close
      fclose(EDStream_);
      break;
  }
}
#endif	/* ED */

void rtFuncEntry_(void) {
  int ii,jj;
  char buf[30];
  if (FuncLevel_<MAXDEP_) {
    SaveStream[FuncLevel_]=MStream_;
    sprintf(buf,"tmp_fun%i",FuncLevel_);
    MStream_=fopen(buf,"w+");
    FuncLevel_++;
    for (jj=0; jj<sizeof(Set); jj++) {
      sf_set[FuncLevel_].CurSet[jj]^=sf_set[FuncLevel_].CurSet[jj];
      uf_set[FuncLevel_].CurSet[jj]^=uf_set[FuncLevel_].CurSet[jj];
      sf_set[FuncLevel_].OrSet[jj]^=sf_set[FuncLevel_].OrSet[jj];
      uf_set[FuncLevel_].OrSet[jj]^=uf_set[FuncLevel_].OrSet[jj];
      for (ii=0; ii<FuncLevel_; ii++) 
	for (jj=0; jj<sizeof(Set); jj++) {
	  sf_set[FuncLevel_].OrSet[jj]|=sf_set[ii].CurSet[jj];
	  uf_set[FuncLevel_].OrSet[jj]|=uf_set[ii].CurSet[jj];
	}
    }
  } else FuncLevel_++;
}

void rtFuncExit_(void) {
  char buf[80];
  FuncLevel_--;
  if (FuncLevel_<MAXDEP_) {
    rewind(MStream_);
    while (fgets(buf,sizeof(buf),MStream_))
      fputs(buf,FStream_);
    fclose(MStream_);
    sprintf(buf,"tmp_fun%i",FuncLevel_);
    remove(buf);
    MStream_=SaveStream[FuncLevel_];
  }
}

static void ResetIndex(void) {
  MemCopy(sf_set[FuncLevel_].CurSet,sf_set[FuncLevel_].OrSet,sizeof(Set));
  MemCopy(uf_set[FuncLevel_].CurSet,uf_set[FuncLevel_].OrSet,sizeof(Set));
}

static int GetIndex(SetOfIndexPtr s) {
  int ii,jj,kk,n;
  unsigned char c;
  for (ii=0; ii<sizeof(s->CurSet); ii++) {
    c=s[FuncLevel_].CurSet[ii] | s[FuncLevel_].OrSet[ii];
    if (c!=0xff) {
      for (jj=0; c&1; jj++,c>>=1);
      s[FuncLevel_].CurSet[ii]|=1<<jj;
      for (kk=0; kk<FuncLevel_; kk++) {
	s[kk].OrSet[ii]|=1<<jj;
	s[kk].CurSet[ii]|=1<<jj;
      }
      n=ii*8+jj;
      if (n+1>s[0].Max) s[0].Max=n+1;
      return n;
    }
  }
  // all are reserved
  return -1;	// let it to crash
}

// Get Index for user-defined functions
int rtGetIndex_(void) {
  return GetIndex(uf_set);
}


// Get Index for standard functions
static int sfGetIndex(void) {
  return GetIndex(sf_set);
}

int rtGetMaxsf_(void) {
  return sf_set[0].Max;
}

int rtGetMaxuf_(void) {
#if DEB
    d_check();
#endif
  return uf_set[0].Max;
}

//
void MemFree__(void *p) {
  MemFree(p);
}

char *Add__(char *op1, char *op2, char freemem) {
  return Add(op1,op2,freemem);
}

char *Mult__(char *op1, char *op2, char freemem) {
  return Mult(op1,op2,freemem);
}