// Enable the External Differentiator
#define ED	1	// MUST be the same as in common.h

typedef enum {DEF,INDEP,TEMP,GRAD} VarKind;	// Kind of variable

class Gradient {
  char **Expr;		// expression(s)
  char ***Der1;		// and gradient(s)
  char ***Der2;		// and second order derivatives
  char ***Der3;		// and third order derivatives
  int n2,n3;		// numbers of 2-d and 3-d order derivatives
  int Vardim;
  int Index;
  VarKind Kind;			// independent, temporary, or result var
  Gradient *Parent;
 public:
  Gradient(char *name=NULL, int varno=0, int vardim=1, VarKind kind=DEF);
  Gradient(double c);
  Gradient(int i);
  Gradient(const Gradient &g);	// copy
  ~Gradient();
  Gradient operator[](int index) {
    Gradient t=*this;
    t.Index=index;
    return t;
  };

  Gradient &operator=(const Gradient &f);
  Gradient &operator=(double c);
  Gradient &operator=(int i);
  Gradient &operator+=(const Gradient &f);
  Gradient &operator+=(double c);
  Gradient &operator+=(int i);
  Gradient &operator-=(const Gradient &f);
  Gradient &operator-=(double c);
  Gradient &operator-=(int i);
  Gradient &operator*=(const Gradient &f);
  Gradient &operator*=(double c);
  Gradient &operator*=(int i);
  Gradient &operator/=(const Gradient &f);
  Gradient &operator/=(double c);
  Gradient &operator/=(int i);

  friend Gradient operator+(const Gradient &f, const Gradient &g);
  friend Gradient operator+(const Gradient &f, double c);
  friend Gradient operator+(double c, const Gradient &g);
  friend Gradient operator+(const Gradient &f, int i);
  friend Gradient operator+(int i, const Gradient &g);

  friend Gradient operator*(const Gradient &f, const Gradient &g);
  friend Gradient operator*(const Gradient &f, double c);
  friend Gradient operator*(double c, const Gradient &g);
  friend Gradient operator*(const Gradient &f, int i);
  friend Gradient operator*(int i, const Gradient &g);

  friend Gradient operator-(const Gradient &f, const Gradient &g);
  friend Gradient operator-(const Gradient &f);
  friend Gradient operator-(const Gradient &f, double c);
  friend Gradient operator-(double c, const Gradient &g);
  friend Gradient operator-(const Gradient &f, int i);
  friend Gradient operator-(int i, const Gradient &g);

  friend Gradient operator/(const Gradient &f, const Gradient &g);
  friend Gradient operator/(const Gradient &f, double c);
  friend Gradient operator/(double c, const Gradient &g);
  friend Gradient operator/(const Gradient &f, int i);
  friend Gradient operator/(int i, const Gradient &g);

  friend Gradient operator>(const Gradient &f, const Gradient &g);
  friend Gradient operator>(const Gradient &f, double c);
  friend Gradient operator>(const Gradient &f, int i);
  friend Gradient operator>(double c, const Gradient &f);
  friend Gradient operator>(int i, const Gradient &f);

  friend Gradient operator>=(const Gradient &f, const Gradient &g);
  friend Gradient operator>=(const Gradient &f, double c);
  friend Gradient operator>=(const Gradient &f, int i);
  friend Gradient operator>=(double c, const Gradient &f);
  friend Gradient operator>=(int i, const Gradient &f);

  friend Gradient operator<(const Gradient &f, const Gradient &g);
  friend Gradient operator<(const Gradient &f, double c);
  friend Gradient operator<(const Gradient &f, int i);
  friend Gradient operator<(double c, const Gradient &f);
  friend Gradient operator<(int i, const Gradient &f);

  friend Gradient operator<=(const Gradient &f, const Gradient &g);
  friend Gradient operator<=(const Gradient &f, double c);
  friend Gradient operator<=(const Gradient &f, int i);
  friend Gradient operator<=(double c, const Gradient &f);
  friend Gradient operator<=(int i, const Gradient &f);

  friend Gradient operator==(const Gradient &f, const Gradient &g);
  friend Gradient operator==(const Gradient &f, double c);
  friend Gradient operator==(const Gradient &f, int i);
  friend Gradient operator==(double c, const Gradient &f);
  friend Gradient operator==(int i, const Gradient &f);

  friend Gradient operator!=(const Gradient &f, const Gradient &g);
  friend Gradient operator!=(const Gradient &f, double c);
  friend Gradient operator!=(const Gradient &f, int i);
  friend Gradient operator!=(double c, const Gradient &f);
  friend Gradient operator!=(int i, const Gradient &f);

  friend Gradient operator!(const Gradient &f);
  friend Gradient operator&&(const Gradient &f, const Gradient &g);
  friend Gradient operator||(const Gradient &f, const Gradient &g);

  friend Gradient fabs(const Gradient &f);
  friend Gradient exp(const Gradient &f);
  friend Gradient log(const Gradient &f);
  friend Gradient sqrt(const Gradient &f);
  friend Gradient pow(const Gradient &f, double p);
  friend Gradient pow(const Gradient &f, int i);
  friend Gradient pow(const Gradient &f, const Gradient &g);
  friend Gradient sin(const Gradient &f);
  friend Gradient cos(const Gradient &f);
  friend Gradient tan(const Gradient &f);
  friend Gradient asin(const Gradient &f);
  friend Gradient acos(const Gradient &f);
  friend Gradient atan(const Gradient &f);
  friend Gradient iff(const Gradient &c, const Gradient &f, const Gradient &g);
#define sqr(x) pow(x,2)
#define cube(x) pow(x,3)
  friend char *GetExpr_(const Gradient &f);
  friend char *GetExprT_(const Gradient &f);
  void SetExpr(char *newtext);
  friend char *GetDer1_(const Gradient &f, int i);
  friend char *GetDer1T_(const Gradient &f, int i);
  void SetDer1(int i, char *newtext);
  friend char *GetDer2_(const Gradient &f, int ic, int jc);
  friend char *GetDer2T_(const Gradient &f, int ic, int jc);
  void SetDer2(int ic, int jc, char *newtext);
  friend char *GetDer3_(const Gradient &f, int ic, int jc, int kc);
  friend char *GetDer3T_(const Gradient &f, int ic, int jc, int kc);
  void SetDer3(int ic, int jc, int kc, char *newtext);
};

#define IndepVar(name,varno,vardim) Gradient name(#name,varno,vardim,INDEP)
#define TempVar(name,vardim)	Gradient name(#name,0,vardim,TEMP)
#define GradVar(name,varno,vardim)  Gradient name(#name,0,vardim,GRAD)

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

#define ind2_(i,j) ((i)+(((j)*((j)+1))>>1))
#define ind3_(i,j,k) ((i)+(((j)*((j)+1))>>1)+(k)*((k)+1)*((k)+2)/6)

#define ind4_(i,j,k,l) (ind3_(i,j,k)+(l)*((l)+1)*((l)+2)*((l)+3)/24)
#define ind5_(i,j,k,l,m) (ind4_(i,j,k,l)+(m)*((m)+1)*((m)+2)*((m)+3)*((m)+4)/120)


extern "C" void rtInit_(FILE *m, FILE *f);
extern "C" void rtTerm_(void);
extern "C" void rtFuncEntry_(void);
extern "C" void rtFuncExit_(void);
extern "C" int rtGetIndex_(void);
extern "C" int rtGetMaxsf_(void);
extern "C" int rtGetMaxuf_(void);

#if ED
extern "C" void rtExtDiff_(char *mio, int PreProc);
#endif	/* ED */

extern "C" void MemFree__(void *p);
extern "C" char *Add__(char *op1, char *op2, char freemem);
extern "C" char *Mult__(char *op1, char *op2, char freemem);