From 520274a73d4192ee41103c0f0843ab2709516da5 Mon Sep 17 00:00:00 2001
From: Xavier Del Campo Romero <xavi.dcr@tutanota.com>
Date: Sun, 14 Apr 2024 21:53:02 +0200
Subject: [PATCH] Import upstream mathexpr

---
 mathexpr/LICENSE                 |   23 +
 mathexpr/mathexpr/example.cpp    |   84 ++
 mathexpr/mathexpr/example_c.cpp  |   86 ++
 mathexpr/mathexpr/mathexpr.cpp   | 1218 ++++++++++++++++++++++++++
 mathexpr/mathexpr/mathexpr.h     |  146 ++++
 mathexpr/mathexpr/mathexpr_c.cpp | 1399 ++++++++++++++++++++++++++++++
 mathexpr/mathexpr/mathexpr_c.h   |  160 ++++
 7 files changed, 3116 insertions(+)
 create mode 100644 mathexpr/LICENSE
 create mode 100644 mathexpr/mathexpr/example.cpp
 create mode 100644 mathexpr/mathexpr/example_c.cpp
 create mode 100644 mathexpr/mathexpr/mathexpr.cpp
 create mode 100644 mathexpr/mathexpr/mathexpr.h
 create mode 100644 mathexpr/mathexpr/mathexpr_c.cpp
 create mode 100644 mathexpr/mathexpr/mathexpr_c.h

diff --git a/mathexpr/LICENSE b/mathexpr/LICENSE
new file mode 100644
index 0000000..f8e0206
--- /dev/null
+++ b/mathexpr/LICENSE
@@ -0,0 +1,23 @@
+Copied from http://www.yann-ollivier.org/mathlib/mathexpr:
+
+This software is available under the MIT Licence:
+
+Copyright (c) 1997-2000 Yann OLLIVIER
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/mathexpr/mathexpr/example.cpp b/mathexpr/mathexpr/example.cpp
new file mode 100644
index 0000000..d6d98f5
--- /dev/null
+++ b/mathexpr/mathexpr/example.cpp
@@ -0,0 +1,84 @@
+/*
+
+example.cpp , example file for mathexpr.cpp 2.0
+
+Copyright (c) 1997-2000 Yann OLLIVIER
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+*/
+
+#include <stdio.h>
+#include "mathexpr.h"
+
+// Compilation :
+// g++ mathexpr.cpp example.cpp -lm
+
+int main()
+{
+
+ // The number that the variable "x" will point to
+ double x;
+ // Creates a variable named "x" and which value will be x
+ RVar xvar ( "x" , &x );
+
+ // Asks for a fomula depending on the variable x, e.g. "sin 2x"
+ char s[500]="";
+ printf("Enter a formula depending on the variable x:\n");
+ gets(s);
+
+ // Creates an operation with this formula. The operation depends on one
+ // variable, which is xvar; the third argument is an array of pointers
+ // to variables; the previous argument is its size
+ RVar* vararray[1]; vararray[0]=&xvar;
+ ROperation op ( s, 1, vararray );
+
+ // Affects (indirectly) a value to xvar
+ x=3;
+ // Printfs the value of the formula for x=3;
+ printf("%s = %G for x=3\n\n", op.Expr(), op.Val() );
+
+ // Creates a function name which can be used in later functions to refer to the operation op(x)
+ RFunction f (op, &xvar); f.SetName("f");
+
+ // Creates a second variable named y, and a formula depending on both x and y
+ double y;
+ RVar yvar ( "y" , &y );
+ RVar* vararray2[2]; // table of variables containing the adresses of xvar and yvar
+ vararray2[0]=&xvar; vararray2[1]=&yvar;
+
+ // Asks for a formula using x, y and the already-defined function f, e.g. x+f(3y)
+ printf("Enter a formula using x, y and the function f(x): x -> %s that you just entered, e.g. x+f(3y) :\n", op.Expr() );
+ gets(s);
+ RFunction* funcarray[1]; funcarray[0]=&f;
+ ROperation op2 ( (char*)s , 2 , vararray2 , 1, funcarray );
+ // vararray2 is a RVar* array with two elements
+ // funcarray is a RFunction* array with one element
+ y=5;
+ printf("Value for x=3, y=5 : %G\n", op2.Val() );
+
+ // Turns the last expression into a function of x and y
+ RFunction g(op2, 2, vararray2); g.SetName("g");
+
+ // Here is another way to do it
+ double z,t;
+ RVar zvar("z", &z), tvar("t", &t);
+ ROperation op3,zop,top;
+ zop=zvar; top=tvar; // constructs, from a variable, the operation returning its value
+
+ op3=g( (zop+top, top^2) ); // Ready-to-use ; needs two pairs of ( )
+ // Now op3 contains the operation op2 with x replaced with z+t, and y replaced with t^2
+
+ z=5;t=7;
+ printf("\nLet g be the function g : x,y -> %s\n", op2.Expr() );
+ printf("Value of %s for z=5,t=7:\n %G\n", op3.Expr(), op3.Val() );
+
+ ROperation dopdt = op3.Diff(tvar); // Computes the derivative of op3 w.r.t t
+ printf("Value of d/dt (g(z+t,t^2)) = %s for z=5,t=7:\n %G\n", dopdt.Expr(), dopdt.Val() );
+
+ return 0;
+}
diff --git a/mathexpr/mathexpr/example_c.cpp b/mathexpr/mathexpr/example_c.cpp
new file mode 100644
index 0000000..f33cf61
--- /dev/null
+++ b/mathexpr/mathexpr/example_c.cpp
@@ -0,0 +1,86 @@
+/*
+
+example.cpp , example file for mathexpr.cpp 2.0
+
+Copyright (c) 1997-2000 Yann OLLIVIER
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+*/
+
+#include <stdio.h>
+#include "mathexpr_c.h"
+
+// Compilation :
+// g++ mathexpr_c.cpp example_c.cpp
+
+int main()
+{
+
+ // The number that the variable "x" will point to
+ double_complex x;
+ // Creates a variable named "x" and which value will be x
+ CVar xvar ( "x" , &x );
+
+ // Asks for a fomula depending on the variable x, e.g. "sin 2x"
+ char s[500]="";
+ printf("Enter a formula depending on the variable x:\n");
+ gets(s);
+
+ // Creates an operation with this formula. The operation depends on one
+ // variable, which is xvar ; the third argument is an array of pointers
+ // to variables; the previous argument is its size
+ CVar* vararray[1]; vararray[0]=&xvar;
+ COperation op ( s, 1, vararray );
+
+ // Affects (indirectly) a value to xvar
+ x=3;
+ // Printfs the value of the formula for x=3;
+ printf("%s = %s for x=3\n\n", op.Expr(), PrettyPrint(op.Val()) );
+
+ // Creates a function name which can be used in later functions to refer to the operation op(x)
+ CFunction f (op, &xvar); f.SetName("f");
+
+ // Creates a second variable named y, and a formula depending on both x and y
+ double_complex y;
+ CVar yvar ( "y" , &y );
+ CVar* vararray2[2]; // table of variables containing the adresses of xvar and yvar
+ vararray2[0]=&xvar; vararray2[1]=&yvar;
+
+ // Asks for a formula using x, y and the already-defined function f, e.g. x+f(3y)
+ printf("Enter a formula using x, y and the function f(x): x -> %s that you just entered, e.g. x+f(3y) :\n", op.Expr());
+ gets(s);
+ CFunction* funcarray[1]; funcarray[0]=&f;
+ COperation op2 ( (char*)s , 2 , vararray2 , 1, funcarray );
+ // vararray2 is a CVar* array with two elements
+ // funcarray is a CFunction* array with one element
+ y=5;
+ printf("Value for x=3, y=5 : %s\n", PrettyPrint(op2.Val()) );
+
+ // Turns the last expression into a function of x and y
+ CFunction g(op2, 2, vararray2); g.SetName("g");
+
+ // Here is another way to do it
+ double_complex z,t;
+ CVar zvar("z", &z), tvar("t", &t);
+ COperation op3,zop,top;
+ zop=zvar; top=tvar; // constructs, from a variable, the operation returning its value
+
+ op3=g( (zop+top, top^2) ); // Ready-to-use ; needs two pairs of ( )
+ // Now op3 contains the operation op2 with x replaced with z+t, and y replaced with t^2
+
+ z=5;t=7;
+ printf("\nLet g be the function g : x,y -> %s\n", op2.Expr());
+ printf("Value of %s for z=5,t=7:\n %s\n", op3.Expr(), PrettyPrint(op3.Val()) );
+
+ COperation dopdt=op3.Diff(tvar); // Computes the derivative of op3 w.r.t t
+ printf("Value of d/dt (g(z+t,t^2)) = %s for z=5,t=7:\n %s\n", dopdt.Expr(), PrettyPrint(dopdt.Val()) );
+ COperation dopdtbar=op3.DiffConj(tvar); // Computes the derivative of op3 w.r.t the conjugate of t
+ printf("Value of d/dtbar (g(z+t,t^2)) = %s for z=5,t=7:\n %s\n", dopdtbar.Expr(), PrettyPrint(dopdtbar.Val()) );
+
+ return 0;
+}
diff --git a/mathexpr/mathexpr/mathexpr.cpp b/mathexpr/mathexpr/mathexpr.cpp
new file mode 100644
index 0000000..885d34e
--- /dev/null
+++ b/mathexpr/mathexpr/mathexpr.cpp
@@ -0,0 +1,1218 @@
+/*
+
+mathexpr.cpp version 2.0
+
+Copyright (c) 1997-2000 Yann OLLIVIER
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+*/
+
+#include"mathexpr.h"
+
+char* MidStr(const char*s,int i1,int i2)
+{
+ if(i1<0||i2>=(int)strlen(s)||i1>i2){
+	 char* cp = new char[1];
+	 cp[0] = '\0';
+	 return cp;
+ }
+ char*s1=new char[i2-i1+2];
+ int i;
+ for(i=i1;i<=i2;i++)s1[i-i1]=s[i];
+ s1[i2-i1+1]=0;return s1;
+}
+
+char* CopyStr(const char*s)
+{char*s1=new char[strlen(s)+1];char*s12=s1;const char*s2=s;
+while((*s12++=*s2++));return s1;}
+
+void InsStr(char*&s,int n,char c)// Warning : deletes the old string
+{if(n<0||n>(int)strlen(s))return;
+char*s1=new char[strlen(s)+2];
+int i;
+for(i=0;i<n;i++)s1[i]=s[i];
+s1[n]=c;for(i=n+1;s[i-1];i++)s1[i]=s[i-1];
+s1[i]=0;
+delete[]s;s=s1;
+}
+
+signed char EqStr(const char*s,const char*s2)
+{if(strlen(s)!=strlen(s2))return 0;
+int i;
+for(i=0;s[i];i++)if(s[i]!=s2[i])return 0;
+return 1;
+}
+
+signed char CompStr(const char*s,int n,const char*s2)
+{if(n<0||n>=(int)strlen(s)||n+(int)strlen(s2)>(int)strlen(s))return 0;
+int i;
+for(i=0;s2[i];i++)if(s[i+n]!=s2[i])return 0;
+return 1;
+}
+
+void DelStr(char*&s,int n)//Deletes the old string
+{char*s1=new char[strlen(s)];
+int i;
+for(i=0;i<n;i++)s1[i]=s[i];
+for(i=n;s[i+1];i++)s1[i]=s[i+1];
+s1[i]=0;
+delete[]s;s=s1;
+}
+
+RVar::RVar(const RVar & rvarp)
+{if(this==&rvarp)return;pval=rvarp.pval;name=CopyStr(rvarp.name);
+}
+
+RVar::RVar(const char*namep,double*pvalp)
+{pval=pvalp;name=CopyStr(namep);}
+
+RVar::~RVar()
+{if(name!=NULL)delete[] name;}
+
+RFunction::RFunction()
+{
+  type=-1;name=new char[1];name[0]=0;
+  nvars=0;ppvar=NULL;pfuncval=NULL;op=ErrVal;buf=NULL;
+}
+
+RFunction::RFunction(double ((*pfuncvalp)(double)))
+{
+  type=0;pfuncval=pfuncvalp;name=new char[1];name[0]=0;
+  nvars=1;ppvar=NULL;op=ErrVal;buf=NULL;
+}
+
+RFunction::RFunction(const RFunction& rfunc)
+{
+  if(this==&rfunc)return;
+  type=rfunc.type;op=rfunc.op;
+  pfuncval=rfunc.pfuncval;
+  name=CopyStr(rfunc.name);
+  nvars=rfunc.nvars;
+  if(rfunc.ppvar!=NULL&&nvars){
+    ppvar=new PRVar[nvars];
+    int i;for(i=0;i<nvars;i++)ppvar[i]=rfunc.ppvar[i];
+    buf=new double[nvars];
+  }else {ppvar=NULL;buf=NULL;}
+}
+
+RFunction::RFunction(const ROperation& opp,RVar* pvarp):op(opp)
+{
+  type=1;name=new char[1];name[0]=0;
+  nvars=1;ppvar=new PRVar[1];ppvar[0]=pvarp;buf=new double[1];
+}
+
+RFunction::RFunction(const ROperation& opp, int nvarsp,RVar**ppvarp):op(opp)
+{
+  type=1;name=new char[1];name[0]=0;
+  nvars=nvarsp;
+  if(nvars){
+    ppvar=new PRVar[nvars];
+    int i;for(i=0;i<nvars;i++)ppvar[i]=ppvarp[i];
+    buf=new double[nvars];
+  }else {ppvar=NULL;buf=NULL;}
+}
+
+RFunction::~RFunction()
+{
+  if(name!=NULL)delete[]name;
+  if(ppvar!=NULL)delete[]ppvar;
+  if(buf!=NULL)delete[]buf;
+}
+
+RFunction& RFunction::operator=(const RFunction& rfunc)
+{
+  if(this==&rfunc)return *this;
+  type=rfunc.type;op=rfunc.op;
+  pfuncval=rfunc.pfuncval;
+  delete[]name;
+  name=CopyStr(rfunc.name);
+  if(ppvar!=NULL)delete[]ppvar;
+  ppvar=NULL;
+  if(buf!=NULL)delete[]buf;
+  buf=NULL;
+  nvars=rfunc.nvars;
+  if(type==1&&nvars){
+    ppvar=new PRVar[nvars];buf=new double[nvars];
+    int i;for(i=0;i<nvars;i++)ppvar[i]=rfunc.ppvar[i];
+  }
+  return *this;
+}
+
+void RFunction::SetName(const char*s)
+{if(name!=NULL)delete[]name;name=CopyStr(s);}
+
+double RFunction::Val(double x) const
+{
+  if(type==-1||nvars>=2)return ErrVal;
+  if(type==0)return (*pfuncval)(x);
+  double xb=*(*ppvar)->pval,y;
+  *(*ppvar)->pval=x;  // Warning : could cause trouble if this value is used in a parallel process
+  y=op.Val();
+  *(*ppvar)->pval=xb;
+  return y;
+}
+
+double RFunction::Val(double*pv) const
+{
+  if(type==-1)return ErrVal;
+  if(type==0)return (*pfuncval)(*pv);
+  double y;
+  int i;
+  for(i=0;i<nvars;i++){
+    buf[i]=*ppvar[i]->pval;
+    // Warning : could cause trouble if this value is used in a parallel process
+    *ppvar[i]->pval=pv[i];
+  }
+  y=op.Val();
+  for(i=0;i<nvars;i++)*ppvar[i]->pval=buf[i];
+  return y;
+}
+
+ROperation::ROperation()
+{op=ErrOp;mmb1=NULL;mmb2=NULL;ValC=ErrVal;pvar=NULL;pvarval=NULL;pfunc=NULL;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;BuildCode();}
+
+ROperation::~ROperation()
+{
+  Destroy();
+}
+
+ROperation::ROperation(const ROperation&ROp)
+{
+  op=ROp.op;pvar=ROp.pvar;pvarval=ROp.pvarval;ValC=ROp.ValC;pfunc=ROp.pfunc;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  if(ROp.mmb1!=NULL)mmb1=new ROperation(*(ROp.mmb1));else mmb1=NULL;
+  if(ROp.mmb2!=NULL)mmb2=new ROperation(*(ROp.mmb2));else mmb2=NULL;
+  BuildCode();
+}
+
+ROperation::ROperation(double x)
+{
+  if(x==ErrVal){op=ErrOp;mmb1=NULL;mmb2=NULL;ValC=ErrVal;}
+  else if(x>=0){op=Num;mmb1=NULL;mmb2=NULL;ValC=x;}
+  else{op=Opp;mmb1=NULL;mmb2=new ROperation(-x);ValC=ErrVal;}
+  pvar=NULL;pvarval=NULL;pfunc=NULL;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  BuildCode();
+}
+
+ROperation::ROperation(const RVar&varp)
+{op=Var;mmb1=NULL;mmb2=NULL;ValC=ErrVal;pvar=&varp;pvarval=varp.pval;containfuncflag=0;pfunc=NULL;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;BuildCode();}
+
+ROperation& ROperation::operator=(const ROperation& ROp)
+{
+  if(this==&ROp)return *this;
+  Destroy();
+  op=ROp.op;pvar=ROp.pvar;pvarval=ROp.pvarval;ValC=ROp.ValC;pfunc=ROp.pfunc;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  if(ROp.mmb1!=NULL)mmb1=new ROperation(*(ROp.mmb1));else mmb1=NULL;
+  if(ROp.mmb2!=NULL)mmb2=new ROperation(*(ROp.mmb2));else mmb2=NULL;
+  BuildCode();
+  return *this;
+}
+
+int operator==(const ROperation& op,const double v)
+{return(op.op==Num&&op.ValC==v);}
+
+int operator==(const ROperation& op1,const ROperation& op2)
+{if(op1.op!=op2.op)return 0;
+if(op1.op==Var)return(*(op1.pvar)==*(op2.pvar));
+if(op1.op==Fun)return(op1.pfunc==op2.pfunc); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
+if(op1.op==Num)return(op1.ValC==op2.ValC);
+if(op1.mmb1==NULL&&op2.mmb1!=NULL)return 0;
+if(op1.mmb2==NULL&&op2.mmb2!=NULL)return 0;
+if(op2.mmb1==NULL&&op1.mmb1!=NULL)return 0;
+if(op2.mmb2==NULL&&op1.mmb2!=NULL)return 0;
+return(((op1.mmb1==NULL&&op2.mmb1==NULL)||(*(op1.mmb1)==*(op2.mmb1)))&&
+	((op1.mmb2==NULL&&op2.mmb2==NULL)||(*(op1.mmb2)==*(op2.mmb2))));
+}
+
+int operator!=(const ROperation& op1,const ROperation& op2)
+{
+  if(op1.op!=op2.op)return 1;
+  if(op1.op==Var)return(op1.pvar!=op2.pvar);
+  if(op1.op==Fun)return(!(op1.pfunc==op2.pfunc)); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
+  if(op1.op==Num)return(op1.ValC!=op2.ValC);
+  if(op1.mmb1==NULL&&op2.mmb1!=NULL)return 1;
+  if(op1.mmb2==NULL&&op2.mmb2!=NULL)return 1;
+  if(op2.mmb1==NULL&&op1.mmb1!=NULL)return 1;
+  if(op2.mmb2==NULL&&op1.mmb2!=NULL)return 1;
+  return(((op1.mmb1!=NULL||op2.mmb1!=NULL)&&(*(op1.mmb1)!=*(op2.mmb1)))||
+	 ((op1.mmb2!=NULL||op2.mmb2!=NULL)&&(*(op1.mmb2)!=*(op2.mmb2))));
+}
+
+ROperation ROperation::operator+() const
+{return *this;}
+
+ROperation ROperation::operator-() const
+{if(op==Num)return -ValC;
+ROperation resultat;
+if(op==Opp)resultat=*mmb2;else{resultat.op=Opp;resultat.mmb2=new ROperation(*this);};
+return resultat;
+}
+
+ROperation operator,(const ROperation& op1,const ROperation& op2)
+{ROperation resultat;
+resultat.op=Juxt;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
+}
+
+ROperation operator+(const ROperation& op1,const ROperation& op2)
+{
+if(op1.op==Num&&op2.op==Num)return op1.ValC+op2.ValC;
+if(op1==0.)return op2;if(op2==0.)return op1;
+if(op1.op==Opp)return op2-*(op1.mmb2);if(op2.op==Opp)return op1-*(op2.mmb2);
+ROperation resultat;
+resultat.op=Add;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
+}
+
+ROperation operator-(const ROperation& op1,const ROperation& op2)
+{
+if(op1.op==Num&&op2.op==Num)return op1.ValC-op2.ValC;
+if(op1==0.)return -op2;if(op2==0.)return op1;
+if(op1.op==Opp)return -(op2+*(op1.mmb2));if(op2.op==Opp)return op1+*(op2.mmb2);
+ROperation resultat;
+resultat.op=Sub;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
+}
+
+ROperation operator*(const ROperation& op1,const ROperation& op2)
+{
+if(op1.op==Num&&op2.op==Num)return op1.ValC*op2.ValC;
+if(op1==0.||op2==0.)return 0.;
+if(op1==1.)return op2;if(op2==1.)return op1;
+if(op1.op==Opp)return -(*(op1.mmb2)*op2);if(op2.op==Opp)return -(op1**(op2.mmb2));
+ROperation resultat;
+resultat.op=Mult;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
+}
+
+ROperation operator/(const ROperation& op1,const ROperation& op2)
+{if(op1.op==Num&&op2.op==Num)return (op2.ValC?op1.ValC/op2.ValC:ErrVal);
+if(op1==0.0)return 0.;if(op2==1.)return op1;if(op2==0.)return ErrVal;
+if(op1.op==Opp)return -(*(op1.mmb2)/op2);if(op2.op==Opp)return -(op1/(*(op2.mmb2)));
+ROperation resultat;
+resultat.op=Div;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
+}
+
+ROperation operator^(const ROperation& op1,const ROperation& op2)
+{if(op1==0.)return 0.;
+if(op2==0.)return 1.;
+if(op2==1.)return op1;
+ROperation resultat;
+resultat.op=Pow;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
+}               
+                
+ROperation sqrt(const ROperation& op)
+{ROperation rop;rop.op=Sqrt;rop.mmb2=new ROperation(op);return rop;}
+ROperation abs(const ROperation& op)
+{ROperation rop;rop.op=Abs;rop.mmb2=new ROperation(op);return rop;}
+ROperation sin(const ROperation& op)
+{ROperation rop;rop.op=Sin;rop.mmb2=new ROperation(op);return rop;}
+ROperation cos(const ROperation& op)
+{ROperation rop;rop.op=Cos;rop.mmb2=new ROperation(op);return rop;}
+ROperation tan(const ROperation& op)
+{ROperation rop;rop.op=Tg;rop.mmb2=new ROperation(op);return rop;}
+ROperation log(const ROperation& op)
+{ROperation rop;rop.op=Ln;rop.mmb2=new ROperation(op);return rop;}
+ROperation exp(const ROperation& op)
+{ROperation rop;rop.op=Exp;rop.mmb2=new ROperation(op);return rop;}
+ROperation acos(const ROperation& op)
+{ROperation rop;rop.op=Acos;rop.mmb2=new ROperation(op);return rop;}
+ROperation asin(const ROperation& op)
+{ROperation rop;rop.op=Asin;rop.mmb2=new ROperation(op);return rop;}
+ROperation atan(const ROperation& op)
+{ROperation rop;rop.op=Atan;rop.mmb2=new ROperation(op);return rop;}
+
+ROperation ApplyOperator(int n,ROperation**pops,ROperation (*func)(const ROperation&,const ROperation&))
+{
+  if(n<=0)return ErrVal;
+  if(n==1)return *pops[0];
+  if(n==2)return (*func)(*pops[0],*pops[1]);
+  return (*func)(*pops[0],ApplyOperator(n-1,pops+1,func));
+}
+
+ROperation RFunction::operator()(const ROperation& op)
+{
+  /* Code to use to replace explcitly instead of using a pointer to
+     if(nvars!=op.NMembers()||type==-1||type==0)return ErrVal;
+     ROperation op2=*pop;int i;
+     RVar**ppvar2=new PRVar[nvars];char s[11]="";
+     for(i=0;i<nvars;i++){
+     sprintf(s,";var%i;",i);
+     ppvar2[i]=new RVar(s,NULL);
+     op2=op2.Substitute(*ppvar[i],(ROperation)*ppvar2[i]);
+     }
+     for(i=0;i<nvars;i++){
+     op2=op2.Substitute(*ppvar2[i],op.NthMember(i+1));
+     delete ppvar2[i];
+     }
+     delete[]ppvar2;
+     return op2;
+  */
+  ROperation op2;op2.op=Fun;op2.pfunc=this;op2.mmb2=new ROperation(op);
+  return op2;
+}
+
+//Auxiliary string functions
+
+void SupprSpaces(char*&s)//Deletes the old string
+{
+ int i;
+ for(i=0;s[i];i++)if(s[i]==' '||s[i]=='\t'||s[i]=='\n')DelStr(s,i--);
+}
+
+signed char IsNumeric(char c)
+{if(c!='0'&&c!='1'&&c!='2'&&c!='3'&&c!='4'
+	&&c!='5'&&c!='6'&&c!='7'&&c!='8'&&c!='9'&&c!='.')return 0;
+return 1;
+}
+
+signed char IsTNumeric(char *s)
+{int i;for(i=0;i<(int)strlen(s);i++)if(!IsNumeric(s[i]))return 0;return 1;
+}
+
+int SearchCorOpenbracket(char*s,int n)  //Searchs the corresponding bracket of an opening bracket
+{if(n>=(int)strlen(s)-1)return -1;
+int i,c=1;
+for(i=n+1;s[i];i++){
+	if(s[i]=='(')c++;else if(s[i]==')')c--;
+	if(!c)return i;
+  };
+return -1;
+}
+
+int SearchCorClosebracket(char*s,int n)  //Searchs the corresponding bracket of a closing bracket
+{if(n<1)return -1;
+int i,c=1;
+for(i=n-1;i>=0;i--){
+	if(s[i]==')')c++;else if(s[i]=='(')c--;
+	if(!c)return i;
+  };
+return -1;
+}
+
+int SearchOperator(char*s,ROperator op)
+{
+  char opc;
+  switch(op){
+  case ErrOp:case Num:case Var:return -1;
+  case Juxt:opc=',';break;
+  case Add:opc='+';break;
+  case Sub:opc='-';break;
+  case Mult:opc='*';break;
+  case Div:opc='/';break;
+  case Pow:opc='^';break;
+  case NthRoot:opc='#';break;
+  case E10:opc='E';break;        
+  default:return -1;
+  };
+  int i;
+  for(i=(int)strlen(s)-1;i>=0;i--){
+    if(s[i]==opc&&(op!=Sub||i&&s[i-1]==')'))return i;
+    if(s[i]==')'){i=SearchCorClosebracket(s,i);if(i==-1)return -1;};
+  };
+  return -1;
+}
+
+void SimplifyStr(char*&s) //Warning : deletes the old string
+{if(!strlen(s))return;
+char*s1=s,*s2=s+strlen(s);signed char ind=0;
+if(s1[0]=='('&&SearchCorOpenbracket(s1,0)==s2-s1-1){
+	s1++;s2--;ind=1;}
+if(s1==s2)
+{
+	delete[]s;
+	s=new char[1]; // ISO C++ forbids initialization in array new
+	s[0]=0;
+	return;
+}
+if(s1[0]==' '){ind=1;while(s1[0]==' '&&s1<s2)s1++;}
+if(s1==s2)
+{
+	delete[]s;
+	s=new char[1]; // ISO C++ forbids initialization in array new
+	s[0]=0;
+	return;
+}
+if(*(s2-1)==' '){ind=1;while(s2>s1&&*(s2-1)==' ')s2--;}
+*s2=0;
+s1=CopyStr(s1);delete[]s;s=s1;
+if(ind)SimplifyStr(s);
+}
+
+int max(int a, int b){return (a>b?a:b);}
+
+int IsVar(const char*s,int n,int nvar,PRVar*ppvar)
+{
+  if(n<0||n>(int)strlen(s))return 0;
+  int i;int l=0;
+  for(i=0;i<nvar;i++)if(CompStr(s,n,(*(ppvar+i))->name))l=max(l,strlen((*(ppvar+i))->name));
+  return l;
+}
+
+int IsFunction(const char*s,int n)
+{
+  if(CompStr(s,n,"sin")||CompStr(s,n,"cos")||CompStr(s,n,"exp")
+     ||CompStr(s,n,"tan")||CompStr(s,n,"log")||CompStr(s,n,"atg")
+     ||CompStr(s,n,"abs"))return 3;
+  if(CompStr(s,n,"tg")||CompStr(s,n,"ln"))return 2;
+  if(CompStr(s,n,"sqrt")||CompStr(s,n,"asin")||CompStr(s,n,"atan")||
+     CompStr(s,n,"acos"))return 4;
+  if(CompStr(s,n,"arcsin")||CompStr(s,n,"arccos")||CompStr(s,n,"arctan"))return 6;
+  if(CompStr(s,n,"arctg"))return 5;
+  return 0;
+}
+
+int IsFunction(const char*s,int n,int nfunc,PRFunction*ppfunc)
+	//Not recognized if a user-defined function is eg "sine" ie begins like
+	//a standard function
+	//IF PATCHED TO DO OTHERWISE, SHOULD BE PATCHED TOGETHER WITH THE
+	//PARSER BELOW which treats standard functions before user-defined ones
+{
+  int l=IsFunction(s,n);
+  if(l)return l;
+  int i;l=0;
+  for(i=0;i<nfunc;i++)if(CompStr(s,n,ppfunc[i]->name))l=max(l,strlen(ppfunc[i]->name));
+  return l;
+}
+
+signed char IsFunction(ROperator op)
+{return (op==Exp||op==Abs||op==Sin||op==Cos||op==Tg||op==Ln||
+	op==Atan||op==Asin||op==Acos||op==Atan||op==Sqrt||op==Opp);
+}
+
+void IsolateVars(char*&s,int nvar,PRVar*ppvar,int nfunc,PRFunction*ppfunc)//Deletes the old string
+{
+  int i,j;
+  i=0;
+  for(i=0;s[i];i++){
+    if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)return;continue;};
+    if(((j=IsVar(s,i,nvar,ppvar))>IsFunction(s,i,nfunc,ppfunc))||((CompStr(s,i,"pi")||CompStr(s,i,"PI")||CompStr(s,i,"Pi"))&&(j=2))){
+      InsStr(s,i,'(');InsStr(s,i+j+1,')');i+=j+1;continue;};
+    if(IsFunction(s,i,nfunc,ppfunc)){i+=IsFunction(s,i,nfunc,ppfunc)-1;if(!s[i])return;continue;};
+  };
+}
+
+void IsolateNumbers(char*&s,int nvar,RVar**ppvar,int nfunc,RFunction**ppfunc)//Deletes the old string
+{
+  int i,i2,ind=0,t1,t2;
+  for(i=0;s[i];i++){
+    if(ind&&!IsNumeric(s[i])){ind=0;InsStr(s,i2,'(');i++;InsStr(s,i,')');continue;};
+		t1=IsVar(s,i,nvar,ppvar);t2=IsFunction(s,i,nfunc,ppfunc);
+    if(t1||t2){i+=max(t1,t2)-1;continue;};
+    if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)return;continue;};
+    if(!ind&&IsNumeric(s[i])){i2=i;ind=1;};
+  };
+if(ind)InsStr(s,i2,'(');i++;InsStr(s,i,')');
+}
+
+ROperation::ROperation(char*sp,int nvar,PRVar*ppvarp,int nfuncp,PRFunction*ppfuncp)
+{
+  ValC=ErrVal;mmb1=NULL;mmb2=NULL;pvar=NULL;op=ErrOp;pvarval=NULL;containfuncflag=0;pfunc=NULL;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  int i,j,k,l;signed char flag=1;
+  char*s=CopyStr(sp),*s1=NULL,*s2=NULL;
+  SimplifyStr(s);if(!s[0]||!strcmp(s,"Error")){goto fin;}
+  while(s[0]==':'||s[0]==';'){
+    s1=CopyStr(s+1);delete[]s;s=s1;s1=NULL;
+    SimplifyStr(s);if(!s[0]||!strcmp(s,"Error")){goto fin;}
+  }
+  if(IsTNumeric(s)){op=Num;ValC=atof(s);mmb1=NULL;mmb2=NULL;goto fin;};
+  if(EqStr(s,"pi")||EqStr(s,"PI")||EqStr(s,"Pi"))
+    {op=Num;ValC=3.141592653589793238462643383279L;mmb1=NULL;mmb2=NULL;goto fin;};
+	if(IsFunction(s,0,nfuncp,ppfuncp)<IsVar(s,0,nvar,ppvarp))
+		for(i=0;i<nvar;i++)if(EqStr(s,(*(ppvarp+i))->name))
+		{pvar=ppvarp[i];pvarval=pvar->pval;op=Var;mmb1=NULL;mmb2=NULL;goto fin;};
+	for(k=0;s[k];k++){
+    if(s[k]=='('){k=SearchCorOpenbracket(s,k);if(k==-1)break;continue;};
+    if((l=IsFunction(s,k,nfuncp,ppfuncp))&&l>=IsVar(s,k,nvar,ppvarp)){
+      i=k+l;while(s[i]==' ')i++;
+      if(s[i]=='('){
+				j=SearchCorOpenbracket(s,i);
+				if(j!=-1){InsStr(s,i,';');k=j+1;}else break;
+      }else if(s[i]!=':'&&s[i]!=';'){InsStr(s,i,':');k=i;}
+    }
+  }
+	IsolateNumbers(s,nvar,ppvarp,nfuncp,ppfuncp);
+	if(nvar)IsolateVars(s,nvar,ppvarp,nfuncp,ppfuncp);
+  SupprSpaces(s);
+  i=SearchOperator(s,Juxt);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Juxt;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Add);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Add;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Sub);
+  if(i!=-1){
+    s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+    op=Sub;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+    mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  if(s[0]=='-'){s2=MidStr(s,1,strlen(s)-1);
+  op=Opp;mmb1=NULL;
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  for(i=0;s[i];i++){
+		if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)break;continue;};
+		if(IsFunction(s,i,nfuncp,ppfuncp)){
+			k=i+IsFunction(s,i,nfuncp,ppfuncp);while(s[k]==' ')k++;
+			if(s[k]==';'){
+				//	s=DelStr(s,k);
+				j=k;while(s[j]!='(')j++;
+				j=SearchCorOpenbracket(s,j);
+				if(j!=-1){InsStr(s,j,')');InsStr(s,i,'(');i=j+2;}
+			}else if(s[k]==':'){
+				//	s=DelStr(s,k);
+				for(j=k;s[j];j++)
+					if(s[j]=='('){j=SearchCorOpenbracket(s,j);break;}
+				if(j==-1)break;
+				for(j++;s[j];j++){
+					if(s[j]=='('){j=SearchCorOpenbracket(s,j);if(j==-1){flag=0;break;};continue;};
+					if(IsFunction(s,j,nfuncp,ppfuncp))break;
+				}
+				if(flag==0){flag=1;break;}
+				while(j>i&&s[j-1]!=')')j--;if(j<=i+1)break;
+				InsStr(s,i,'(');InsStr(s,j+1,')');
+				i=j+1;
+			}
+		}
+  }
+  for(i=0;s[i]&&s[i+1];i++)if(s[i]==')'&&s[i+1]=='(')
+    InsStr(s,++i,'*');
+  if(s[0]=='('&&SearchCorOpenbracket(s,0)==(int)strlen(s)-1){
+    if(CompStr(s,1,"exp")){op=Exp;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"abs")){op=Abs;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"sin")){op=Sin;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"cos")){op=Cos;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"tan")){op=Tg;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"log")){op=Ln;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"atg")){op=Atan;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"tg")){op=Tg;s2=MidStr(s,3,strlen(s)-2);}
+    else if(CompStr(s,1,"ln")){op=Ln;s2=MidStr(s,3,strlen(s)-2);}
+    else if(CompStr(s,1,"asin")){op=Asin;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"acos")){op=Acos;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"atan")){op=Atan;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"sqrt")){op=Sqrt;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"arcsin")){op=Asin;s2=MidStr(s,7,strlen(s)-2);}
+    else if(CompStr(s,1,"arccos")){op=Acos;s2=MidStr(s,7,strlen(s)-2);}
+    else if(CompStr(s,1,"arctan")){op=Atan;s2=MidStr(s,7,strlen(s)-2);}
+    else if(CompStr(s,1,"arctg")){op=Atan;s2=MidStr(s,6,strlen(s)-2);}
+    else {
+			for(i=-1,k=0,j=0;j<nfuncp;j++)if(CompStr(s,1,ppfuncp[j]->name)&&k<(int)strlen(ppfuncp[j]->name)){k=strlen(ppfuncp[j]->name);i=j;}
+			if(i>-1){
+				op=Fun;s2=MidStr(s,strlen(ppfuncp[i]->name)+1,strlen(s)-2);
+				pfunc=ppfuncp[i];
+			}
+		}
+		mmb1=NULL;mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);
+		if(op==Fun)if(mmb2->NMembers()!=pfunc->nvars){op=ErrOp;mmb1=NULL;mmb2=NULL;goto fin;}
+    goto fin;
+  };
+  i=SearchOperator(s,Mult);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Mult;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Div);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Div;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Pow);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Pow;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,NthRoot);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  if(i==0||s[i-1]!=')')
+    {op=Sqrt;mmb1=NULL;}else
+      {op=NthRoot;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);};
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,E10);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=E10;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  op=ErrOp;mmb1=NULL;mmb2=NULL;
+fin:
+	BuildCode();
+	delete[]s;if(s1!=NULL)delete[] s1;if(s2!=NULL)delete[]s2;
+}
+
+void ROperation::Destroy()
+{
+  if(mmb1!=NULL&&mmb2!=NULL&&mmb1!=mmb2){delete mmb1;delete mmb2;mmb1=NULL;mmb2=NULL;}
+  else if(mmb1!=NULL){delete mmb1;mmb1=NULL;}else if(mmb2!=NULL){delete mmb2;mmb2=NULL;}
+  if(pinstr!=NULL){delete[]pinstr;pinstr=NULL;}
+  if(pvals!=NULL){
+		if(op==ErrOp||op==Num)delete pvals[0];
+		delete[]pvals;pvals=NULL;
+	}
+  if(ppile!=NULL){delete[]ppile;ppile=NULL;}
+  if(pfuncpile!=NULL){delete[]pfuncpile;pfuncpile=NULL;}
+}
+
+int operator==(const RVar& var1,const RVar& var2)
+{return(var1.pval==var2.pval&&EqStr(var1.name,var2.name));
+}
+
+int operator==(const RFunction& f1,const RFunction& f2)
+{
+  if(f1.type!=f2.type)return 0;
+  if(f1.type==-1)return 1; // Nonfunction==nonfunction
+  if(f1.type==0)return (f1.pfuncval==f2.pfuncval&&EqStr(f1.name,f2.name));
+  if(f1.op!=f2.op)return 0;
+  if(!EqStr(f1.name,f2.name))return 0;
+  if(f1.nvars!=f2.nvars)return 0;
+  int i;
+  for(i=0;i<f1.nvars;i++)if(!(*f1.ppvar[i]==*f2.ppvar[i]))return 0;
+  return 1;
+}
+
+/*
+double ROperation::Val() const // Won't work if multi-variable functions are included
+{
+  double v1=ErrVal,v2=ErrVal;
+  if(mmb1!=NULL){v1=mmb1->Val();if(fabsl(v1)<sqrtminfloat)v1=0;else if(v1==ErrVal||fabsl(v1)>sqrtmaxfloat)return ErrVal;};
+  if(mmb2!=NULL){v2=mmb2->Val();if(fabsl(v2)<sqrtminfloat)v2=0;else if(v2==ErrVal||fabsl(v2)>sqrtmaxfloat)return ErrVal;};
+  switch(op){
+  case Num:return ValC;
+  case Var:return *pvarval;
+  case Add:return v1+v2;
+  case Sub:return v1-v2;
+  case Opp:return -v2;
+  case Mult:return v1*v2;
+  case Div:if(v2)return v1/v2;else return ErrVal;
+  case Pow:if(v1==0)return 0;else if((v1>0||!fmodl(v2,1))&&v2*logl(fabsl(v1))<DBL_MAX_EXP)return powl(v1,v2);else return ErrVal;
+  case Sqrt:if(v2>=0)return sqrtl(v2);else return ErrVal;
+  case NthRoot:if(!v1||v2*logl(fabsl(v1))<DBL_MIN_EXP)return ErrVal;
+  else if(v2>=0)return powl(v2,1/v1);else
+    if(fmodl(v1,2)==1||fmodl(v1,2)==-1)return -powl(-v2,1/v1);else return ErrVal;
+  case E10:if(v2<DBL_MAX_10_EXP)return v1*pow10l(v2);else return ErrVal;
+  case Ln:if(v2>0)return logl(v2);else return ErrVal;
+  case Exp:if(v2<DBL_MAX_EXP)return expl(v2);else return ErrVal;
+  case Sin:if(fabsl(v2)<inveps)return sinl(v2);else return ErrVal;
+  case Cos:if(fabsl(v2)<inveps)return cosl(v2);else return ErrVal;
+  case Tg:if(fabsl(v2)<inveps)return tanl(v2);else return ErrVal;
+  case Atan:
+    if(mmb2->op==Juxt){v1=mmb2->NthMember(1).Val();v2=mmb2->NthMember(2).Val();return (v1||v2?atan2(v1,v2):ErrVal);}else return atanl(v2);
+  case Asin:if(v2<-1||v2>1)return ErrVal;else return asinl(v2);
+  case Acos:if(v2<-1||v2>1)return ErrVal;else return acosl(v2);
+  case Abs:return fabsl(v2);
+  case Fun:return pfunc->Val(v2);
+  default:return ErrVal;
+  };
+}
+*/
+
+signed char ROperation::ContainVar(const RVar& varp) const
+{if(op==Var){if(EqStr(pvar->name,varp.name)&&pvar->pval==varp.pval)
+	return 1;else return 0;};
+if(mmb1!=NULL&&mmb1->ContainVar(varp))return 1;
+if(mmb2!=NULL&&mmb2->ContainVar(varp))return 1;
+return 0;
+}
+
+signed char ROperation::ContainFuncNoRec(const RFunction& func) const // No recursive test on subfunctions
+{if(op==Fun){if(*pfunc==func)
+  return 1;else return 0;}
+ if(mmb1!=NULL&&mmb1->ContainFuncNoRec(func))return 1;
+ if(mmb2!=NULL&&mmb2->ContainFuncNoRec(func))return 1;
+ return 0;
+}
+
+signed char ROperation::ContainFunc(const RFunction& func) const // Recursive test on subfunctions
+{
+  if(containfuncflag)return 0;
+  if(op==Fun&&*pfunc==func)return 1;
+  containfuncflag=1;
+  if(op==Fun)if(pfunc->op.ContainFunc(func)){containfuncflag=0;return 1;}
+  if(mmb1!=NULL&&mmb1->ContainFunc(func)){containfuncflag=0;return 1;}
+  if(mmb2!=NULL&&mmb2->ContainFunc(func)){containfuncflag=0;return 1;}
+  containfuncflag=0;return 0;
+}
+
+signed char ROperation::HasError(const ROperation*pop) const
+{
+  if(op==ErrOp)return 1;
+  if(op==Fun&&pfunc->type==1&&pfunc->op==*(pop==NULL?this:pop))return 1;
+  if(op==Fun&&pfunc->type==1&&pfunc->op.HasError((pop==NULL?this:pop)))return 1;
+  if(mmb1!=NULL&&mmb1->HasError((pop==NULL?this:pop)))return 1;
+  if(mmb2!=NULL&&mmb2->HasError((pop==NULL?this:pop)))return 1;
+  if(op==Fun&&pfunc->type==-1)return 1;
+  return 0;
+}
+
+int ROperation::NMembers() const //Number of members for an operation like a,b,c...
+{
+  if(op==Fun)return(pfunc->type==1?pfunc->op.NMembers():pfunc->type==0?1:0);
+  if(op!=Juxt)return 1;else if(mmb2==NULL)return 0;else return 1+mmb2->NMembers();
+}
+
+ROperation ROperation::NthMember(int n) const
+{
+  PRFunction prf;
+  if(op==Fun&&pfunc->type==1&&pfunc->op.NMembers()>1){
+    prf=new RFunction(pfunc->op.NthMember(n),pfunc->nvars,pfunc->ppvar);
+    char*s=new char[strlen(pfunc->name)+10];
+    sprintf(s,"(%s_%i)",pfunc->name,n);prf->SetName(s);delete[]s;
+    return(*prf)(*mmb2);
+  }
+  if(n==1){
+    if(op!=Juxt)return *this; else if(mmb1!=NULL)return *mmb1;else return ErrVal;
+  };
+  if(op!=Juxt)return ErrVal;
+  if(n>1&&mmb2!=NULL)return mmb2->NthMember(n-1);
+  return ErrVal;
+}
+
+ROperation ROperation::Substitute(const RVar& var,const ROperation& rop) const // Replaces variable var with expression rop
+{
+  if(!ContainVar(var))return *this;
+  if(op==Var)return rop;
+  ROperation r;
+  r.op=op;r.pvar=pvar;r.pvarval=pvarval;r.ValC=ValC;r.pfunc=pfunc;
+  if(mmb1!=NULL)r.mmb1=new ROperation(mmb1->Substitute(var,rop));else r.mmb1=NULL;
+  if(mmb2!=NULL)r.mmb2=new ROperation(mmb2->Substitute(var,rop));else r.mmb2=NULL;
+  return r;
+}
+
+ROperation ROperation::Diff(const RVar& var) const
+{
+  if(!ContainVar(var))return 0.0;
+  if(op==Var)return 1.0;
+  ROperation **ppop1,op2;int i,j;
+  switch(op){
+  case Juxt:return(mmb1->Diff(var),mmb2->Diff(var));
+  case Add:return(mmb1->Diff(var)+mmb2->Diff(var));
+  case Sub:return(mmb1->Diff(var)-mmb2->Diff(var));
+  case Opp:return(-mmb2->Diff(var));
+  case Mult:return((*mmb1)*(mmb2->Diff(var))+(*mmb2)*(mmb1->Diff(var)));
+  case Div:if(mmb2->ContainVar(var))return(((*mmb2)*(mmb1->Diff(var))-(*mmb1)*(mmb2->Diff(var)))/((*mmb2)^2));
+  else return(mmb1->Diff(var)/(*mmb2));
+  case Pow:if(mmb2->ContainVar(var))return((*this)*(log(*mmb1)*mmb2->Diff(var)+
+						    (*mmb2)*mmb1->Diff(var)/(*mmb1)));else
+						      return (*mmb2)*mmb1->Diff(var)*((*mmb1)^(*mmb2-1));
+  case Sqrt:return(mmb2->Diff(var)/(2*sqrt(*mmb2)));
+  case NthRoot:{ROperation interm=(*mmb2)^(1/(*mmb1));return interm.Diff(var);};
+  case E10:{ROperation interm=(*mmb1)*(10^(*mmb2));return interm.Diff(var);};;
+  case Ln:return (mmb2->Diff(var)/(*mmb2));
+  case Exp:return (mmb2->Diff(var)*(*this));
+  case Sin:return (mmb2->Diff(var)*cos(*mmb2));
+  case Cos:return (-mmb2->Diff(var)*sin(*mmb2));
+  case Tg:return (mmb2->Diff(var)*(1+((*this)^2)));
+  case Atan:
+    if(mmb2->op!=Juxt)return(mmb2->Diff(var)/(1+((*mmb2)^2)));
+    else return ((mmb2->NthMember(1).Diff(var))*(mmb2->NthMember(2))-(mmb2->NthMember(2).Diff(var))*(mmb2->NthMember(1)))/(((mmb2->NthMember(1))^2)+((mmb2->NthMember(2))^2));
+  case Asin:return(mmb2->Diff(var)/sqrt(1-((*mmb2)^2)));
+  case Acos:return(-mmb2->Diff(var)/sqrt(1-((*mmb2)^2)));
+  case Abs:return(mmb2->Diff(var)*(*mmb2)/(*this));
+  case Fun:if(pfunc->type==-1||pfunc->type==0)return ErrVal;
+    if(pfunc->nvars==0)return 0.;
+    else if(pfunc->op.NMembers()>1){
+      j=pfunc->op.NMembers();
+      ppop1=new ROperation*[j];
+      for(i=0;i<j;i++)ppop1[i]=new ROperation(NthMember(i+1).Diff(var));
+      op2=ApplyOperator(pfunc->nvars,ppop1,&operator,);
+      for(i=0;i<pfunc->nvars;i++)delete ppop1[i];
+      delete[]ppop1;
+      return op2;  
+    }else{
+      ppop1=new ROperation*[pfunc->nvars];
+      for(i=0;i<pfunc->nvars;i++){
+	ppop1[i]=new ROperation(pfunc->op.Diff(*pfunc->ppvar[i]));
+	for(j=0;j<pfunc->nvars;j++)
+	  *ppop1[i]=ppop1[i]->Substitute(*pfunc->ppvar[j],mmb2->NthMember(j+1));
+	*ppop1[i]=(mmb2->NthMember(i+1).Diff(var))*(*ppop1[i]);
+      }
+      op2=ApplyOperator(pfunc->nvars,ppop1,&::operator+);
+      for(i=0;i<pfunc->nvars;i++)delete ppop1[i];
+      delete[]ppop1;
+      return op2;  
+      // In the obtained expression, f' will have been replaced with its expression but f will remain pointing to itself ; this could cause some trouble if changing f afterwards
+    }
+  default:return ErrVal;
+  };
+}
+
+char* ValToStr(double x)
+{
+  char*s=new char[30];
+  if(x==(double)3.141592653589793238462643383279L)sprintf(s,"pi");else sprintf(s,"%.16G",x);
+  return s;
+}
+
+char* ROperation::Expr() const
+{
+  char*s=NULL,*s1=NULL,*s2=NULL;int n=10;signed char f=0,g=0;
+  if(op==Fun)if(strlen(pfunc->name)>4)n+=strlen(pfunc->name)-4;
+  if(mmb1!=NULL){s1=mmb1->Expr();n+=strlen(s1);f=IsFunction(mmb1->op);}
+  if(mmb2!=NULL){s2=mmb2->Expr();n+=strlen(s2);g=IsFunction(mmb2->op);}
+  s=new char[n];
+  switch(op){
+  case Num:return ValToStr(ValC);
+  case Var:return CopyStr(pvar->name);
+  case Juxt:sprintf(s,"%s , %s",s1,s2);break;
+  case Add:
+    f=f||(mmb1->op==Juxt);
+    g=g||(mmb2->op==Juxt);
+    if(f&&g)sprintf(s,"(%s)+(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)+%s",s1,s2);else
+	if(g)sprintf(s,"%s+(%s)",s1,s2);else
+	  sprintf(s,"%s+%s",s1,s2);
+    break;
+  case Sub:
+    f=f||(mmb1->op==Juxt);
+    g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub);
+    if(f&&g)sprintf(s,"(%s)-(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)-%s",s1,s2);else
+	if(g)sprintf(s,"%s-(%s)",s1,s2);else
+	  sprintf(s,"%s-%s",s1,s2);
+    break;
+  case Opp:
+    if(mmb2->op==Add||mmb2->op==Sub||mmb2->op==Juxt)sprintf(s,"-(%s)",s2);else
+      sprintf(s,"-%s",s2);
+    break;
+  case Mult:
+    f=f||(mmb1->op==Juxt||mmb1->op==Add||mmb1->op==Sub||mmb1->op==Opp||mmb1->op==Div);
+    g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub||mmb2->op==Opp);
+    if(f&&g)sprintf(s,"(%s)*(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)*%s",s1,s2);else
+	if(g)sprintf(s,"%s*(%s)",s1,s2);else
+	  sprintf(s,"%s*%s",s1,s2);
+    break;
+  case Div:
+    f=f||(mmb1->op==Juxt||mmb1->op==Add||mmb1->op==Sub||mmb1->op==Opp||mmb1->op==Div);
+    g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub||mmb2->op==Opp||mmb2->op==Mult||mmb2->op==Div);
+    if(f&&g)sprintf(s,"(%s)/(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)/%s",s1,s2);else
+	if(g)sprintf(s,"%s/(%s)",s1,s2);else
+	  sprintf(s,"%s/%s",s1,s2);
+    break;
+  case Pow:
+    f=(mmb1->op!=Num&&mmb1->op!=Var);
+    g=(mmb2->op!=Num&&mmb2->op!=Var);
+    if(f&&g)sprintf(s,"(%s)^(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)^%s",s1,s2);else
+	if(g)sprintf(s,"%s^(%s)",s1,s2);else
+	  sprintf(s,"%s^%s",s1,s2);
+    break;
+  case Sqrt:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"sqrt(%s)",s2);
+    else sprintf(s,"sqrt %s",s2);
+    break;
+  case NthRoot:
+    f=(mmb1->op!=Num&&mmb1->op!=Var);
+    g=(mmb2->op!=Num&&mmb2->op!=Var);
+    if(f&&g)sprintf(s,"(%s)#(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)#%s",s1,s2);else
+	if(g)sprintf(s,"%s#(%s)",s1,s2);else
+	  sprintf(s,"%s#%s",s1,s2);
+    break;
+  case E10:
+    f=(mmb1->op!=Num&&mmb1->op!=Var);
+    g=(mmb2->op!=Num&&mmb2->op!=Var);
+    if(f&&g)sprintf(s,"(%s)E(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)E%s",s1,s2);else
+	if(g)sprintf(s,"%sE(%s)",s1,s2);else
+	  sprintf(s,"%sE%s",s1,s2);
+    break;
+  case Ln:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"log(%s)",s2);
+    else sprintf(s,"log %s",s2);
+    break;
+  case Exp:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"exp(%s)",s2);
+    else sprintf(s,"exp %s",s2);
+    break;
+  case Sin:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"sin(%s)",s2);
+    else sprintf(s,"sin %s",s2);
+    break;
+  case Cos:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"cos(%s)",s2);
+    else sprintf(s,"cos %s",s2);
+    break;
+  case Tg:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"tan(%s)",s2);
+    else sprintf(s,"tan %s",s2);
+    break;
+  case Atan:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"atan(%s)",s2);
+    else sprintf(s,"atan %s",s2);
+    break;
+  case Asin:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"asin(%s)",s2);
+    else sprintf(s,"asin %s",s2);
+    break;
+  case Acos:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"acos(%s)",s2);
+    else sprintf(s,"acos %s",s2);
+    break;
+  case Abs:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"abs(%s)",s2);
+    else sprintf(s,"abs %s",s2);
+    break;
+  case Fun:
+    sprintf(s,"%s(%s)",pfunc->name,s2);
+    break;
+  default:return CopyStr("Error");
+  };
+  if(s1!=NULL)delete[] s1;if(s2!=NULL)delete[] s2;
+  return s;
+}
+
+const double sqrtmaxfloat=sqrt(DBL_MAX);
+const double sqrtminfloat=sqrt(DBL_MIN);
+const double inveps=.1/DBL_EPSILON;
+
+void  Addition(double*&p)
+{if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||fabsl(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p+=(*(p+1));}
+void  Soustraction(double*&p)
+{if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||fabsl(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p-=(*(p+1));}
+void  Multiplication(double*&p)
+{if(fabsl(*p)<sqrtminfloat){*--p=0;return;};
+if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(fabsl(*(--p))<sqrtminfloat){*p=0;return;};
+if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p*=(*(p+1));}
+void  Division(double*&p)
+{if(fabsl(*p)<sqrtminfloat||*p==ErrVal||fabsl(*p)>sqrtmaxfloat)
+	{*(--p)=ErrVal;return;};
+if(fabsl(*(--p))<sqrtminfloat)*p=0;else if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat)
+	{*p=ErrVal;return;};
+*p/=(*(p+1));}
+void  Puissance(double*&p)
+{double v2=*p--,v1=*p;
+if(!v1){*p=0;return;};
+if(v2==ErrVal||v1==ErrVal||fabsl(v2*logl(fabsl(v1)))>DBL_MAX_EXP){*p=ErrVal;return;};
+*p=((v1>0||!fmodl(v2,1))?powl(v1,v2):ErrVal);}
+void  RacineN(double*&p)
+{double v2=*p--,v1=*p;
+if(v1==ErrVal||v2==ErrVal||!v1||v2*logl(fabsl(v1))<DBL_MIN_EXP){*p=ErrVal;return;};
+if(v2>=0){*p=powl(v2,1/v1);return;};
+*p=((fabsl(fmodl(v1,2))==1)?-powl(-v2,1/v1):ErrVal);}
+void  Puiss10(double*&p)
+{if(fabsl(*p)<sqrtminfloat){*(--p)=0;return;};
+if(*p==ErrVal||fabsl(*p)>DBL_MAX_10_EXP){*(--p)=ErrVal;return;};
+if(fabsl(*(--p))<sqrtminfloat)*p=0;else if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat)
+{*p=ErrVal;return;};
+*p*=pow10l(*(p+1));}
+void  ArcTangente2(double*&p)
+{if(*p==ErrVal||fabsl(*p)>inveps){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||fabsl(*p)>inveps){*p=ErrVal;return;};
+*p=(*p||*(p+1)?atan2(*p,*(p+1)):ErrVal);}
+void  NextVal(double*&){}
+void  RFunc(double*&){}
+void  JuxtF(double*&){}
+void  Absolu(double*&p){*p=((*p==ErrVal)?ErrVal:fabsl(*p));}
+void  Oppose(double*&p){*p=((*p==ErrVal)?ErrVal:-*p);}
+void  ArcSinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>1)?ErrVal:asinl(*p));}
+void  ArcCosinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>1)?ErrVal:acosl(*p));}
+void  ArcTangente(double*&p)
+{*p=((*p==ErrVal)?ErrVal:atanl(*p));}
+void  Logarithme(double*&p)
+{*p=((*p==ErrVal||*p<=0)?ErrVal:logl(*p));}
+void  Exponentielle(double*&p)
+{*p=((*p==ErrVal||*p>DBL_MAX_EXP)?ErrVal:expl(*p));}
+void  Sinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>inveps)?ErrVal:sinl(*p));}
+void  Tangente(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>inveps)?ErrVal:tanl(*p));}
+void  Cosinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>inveps)?ErrVal:cosl(*p));}
+void  Racine(double*&p)
+{*p=((*p==ErrVal||*p>sqrtmaxfloat||*p<0)?ErrVal:sqrtl(*p));}
+void FonctionError(double*&p){*p=ErrVal;}
+inline void ApplyRFunc(PRFunction rf,double*&p)
+{p-=rf->nvars-1;*p=rf->Val(p);}
+
+double ROperation::Val() const
+{
+  pfoncld*p1=pinstr;double**p2=pvals,*p3=ppile-1;PRFunction*p4=pfuncpile;
+  for(;*p1!=NULL;p1++)
+    if(*p1==&NextVal)*(++p3)=**(p2++);else 
+      if(*p1==&RFunc) ApplyRFunc(*(p4++),p3);
+      else (**p1)(p3);
+  return *p3;
+}
+
+void BCDouble(pfoncld*&pf,pfoncld*pf1,pfoncld*pf2,
+	      double**&pv,double**pv1,double**pv2,
+	      double*&pp,double*pp1,double*pp2,
+	      RFunction**&prf,RFunction**prf1,RFunction**prf2,
+	      pfoncld f)
+{
+  pfoncld*pf3,*pf4=pf1;long n1,n2;
+  for(n1=0;*pf4!=NULL;pf4++,n1++);for(n2=0,pf4=pf2;*pf4!=NULL;pf4++,n2++);
+  pf=new pfoncld[n1+n2+2];
+  for(pf3=pf,pf4=pf1;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  for(pf4=pf2;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  *pf3++=f;*pf3=NULL;//delete[]pf1,pf2;
+  double**pv3,**pv4=pv1;
+  for(n1=0;*pv4!=NULL;pv4++,n1++);for(n2=0,pv4=pv2;*pv4!=NULL;pv4++,n2++);
+  pv=new double*[n1+n2+1];
+  for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  for(pv4=pv2;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  *pv3=NULL;//delete[]pv1,pv2;
+  double*pp3,*pp4=pp1;
+  for(n1=0;*pp4!=ErrVal;pp4++,n1++);for(n2=0,pp4=pp2;*pp4!=ErrVal;pp4++,n2++);
+  pp=new double[n1+n2+1];  // Really need to add and not to take max(n1,n2) in case of Juxt operator
+  for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  for(pp4=pp2;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  *pp3=ErrVal;//delete[]pp1,pp2;
+  PRFunction*prf3,*prf4=prf1;
+  for(n1=0;*prf4!=NULL;prf4++,n1++);for(n2=0,prf4=prf2;*prf4!=NULL;prf4++,n2++);
+  prf=new PRFunction[n1+n2+1];
+  for(prf3=prf,prf4=prf1;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  for(prf4=prf2;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  *prf3=NULL;//delete[]prf1,prf2;
+}
+
+void BCSimple(pfoncld*&pf,pfoncld*pf1,double**&pv,double**pv1,
+	double*&pp,double*pp1,RFunction**&prf,RFunction**prf1,pfoncld f)
+{
+  pfoncld*pf3,*pf4=pf1;long n;
+  for(n=0;*pf4!=NULL;pf4++,n++);
+  pf=new pfoncld[n+2];
+  for(pf4=pf1,pf3=pf;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  *pf3++=f;*pf3=NULL;//delete[]pf1;
+  double**pv3,**pv4=pv1;
+  for(n=0;*pv4!=NULL;pv4++,n++);
+  pv=new double*[n+1];
+  for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  *pv3=NULL;//delete[]pv1;
+  double*pp3,*pp4=pp1;
+  for(n=0;*pp4!=ErrVal;pp4++,n++);
+  pp=new double[n+1];
+  for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  *pp3=ErrVal;//delete[]pp1;
+  RFunction**prf3,**prf4=prf1;
+  for(n=0;*prf4!=NULL;prf4++,n++);
+  prf=new RFunction*[n+1];
+  for(prf3=prf,prf4=prf1;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  *prf3=NULL;//delete[]prf1;
+}
+
+void BCFun(pfoncld*&pf,pfoncld*pf1,double**&pv,double**pv1,
+	double*&pp,double*pp1,RFunction**&prf,RFunction**prf1,PRFunction rf)
+{
+  pfoncld*pf3,*pf4=pf1;long n;
+  for(n=0;*pf4!=NULL;pf4++,n++);
+  pf=new pfoncld[n+2];
+  for(pf4=pf1,pf3=pf;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  *pf3++=&RFunc;*pf3=NULL;//delete[]pf1;
+  double**pv3,**pv4=pv1;
+  for(n=0;*pv4!=NULL;pv4++,n++);
+  pv=new double*[n+1];
+  for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  *pv3=NULL;//delete[]pv1;
+  double*pp3,*pp4=pp1;
+  for(n=0;*pp4!=ErrVal;pp4++,n++);
+  pp=new double[n+1];
+  for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  *pp3=ErrVal;//delete[]pp1;
+  PRFunction*prf3,*prf4=prf1;
+  for(n=0;*prf4!=NULL;prf4++,n++);
+  prf=new PRFunction[n+2];
+  for(prf4=prf1,prf3=prf;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  *prf3++=rf;*prf3=NULL;//delete[]pf1;
+}
+
+void ROperation::BuildCode()
+{
+  //  if(mmb1!=NULL)mmb1->BuildCode();if(mmb2!=NULL)mmb2->BuildCode();
+  if(pinstr!=NULL){delete[]pinstr;pinstr=NULL;}
+  if(pvals!=NULL){delete[]pvals;pvals=NULL;}//does not delete pvals[0] in case it was to be deleted... (no way to know)
+  if(ppile!=NULL){delete[]ppile;ppile=NULL;}
+  if(pfuncpile!=NULL){delete[]pfuncpile;pfuncpile=NULL;}
+  switch(op){
+  case ErrOp:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+    pvals=new double*[2];pvals[0]=new double(ErrVal);pvals[1]=NULL;
+    ppile=new double[2];ppile[0]=0;ppile[1]=ErrVal;
+    pfuncpile=new RFunction*[1];pfuncpile[0]=NULL;
+    break;
+  case Num:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+    pvals=new double*[2];pvals[0]=new double(ValC);pvals[1]=NULL;
+    ppile=new double[2];ppile[0]=0;ppile[1]=ErrVal;
+    pfuncpile=new RFunction*[1];pfuncpile[0]=NULL;break;
+  case Var:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+    pvals=new double*[2];pvals[0]=pvarval;pvals[1]=NULL;
+    ppile=new double[2];ppile[0]=0;ppile[1]=ErrVal;
+    pfuncpile=new RFunction*[1];pfuncpile[0]=NULL;break;
+  case Juxt:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+		     pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&JuxtF);
+  break;case Add:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+		    pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Addition);
+  break;case Sub:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Soustraction);
+  break;case Mult:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			   pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Multiplication);
+  break;case Div:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Division);
+  break;case Pow:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Puissance);
+  break;case NthRoot:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			      pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&RacineN);
+  break;case E10:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Puiss10);
+  break;case Opp:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Oppose);
+  break;case Sin:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Sinus);
+  break;case Sqrt:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Racine);
+  break;case Ln:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			 ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Logarithme);
+  break;case Exp:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Exponentielle);
+  break;case Cos:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Cosinus);
+  break;case Tg:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			 ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Tangente);
+  break;case Atan:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,(mmb2->NMembers()>1?&ArcTangente2:&ArcTangente));
+  break;case Asin:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&ArcSinus);
+  break;case Acos:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&ArcCosinus);
+  break;case Abs:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Absolu);
+  break;case Fun:BCFun(pinstr,mmb2->pinstr,pvals,mmb2->pvals,ppile,
+		       mmb2->ppile,pfuncpile,mmb2->pfuncpile,pfunc);
+  break;default:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			 ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&FonctionError);
+  }
+}
diff --git a/mathexpr/mathexpr/mathexpr.h b/mathexpr/mathexpr/mathexpr.h
new file mode 100644
index 0000000..4ea9c40
--- /dev/null
+++ b/mathexpr/mathexpr/mathexpr.h
@@ -0,0 +1,146 @@
+/*
+
+mathexpr.h version 2.0
+
+Copyright (c) 1997-2000 Yann OLLIVIER
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+*/
+
+#ifndef _MATHEXPR_H
+#define _MATHEXPR_H
+
+#include<string.h>
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include<float.h>
+
+// Compatibility with long double-typed functions
+#define atanl atan
+#define asinl asin
+#define acosl acos
+#define expl exp
+#define logl log
+#define powl pow
+#define pow10l(x) pow(10,x)
+#define fabsl fabs
+#define cosl cos
+#define sinl sin
+#define tanl tan
+#define fmodl fmod
+#define sqrtl sqrt
+
+
+// Warning : if ints are short, everything will fail with strings longer than 32767 chars
+
+
+const double ErrVal=DBL_MAX;
+
+//Class definitions for operations
+
+class RVar{
+ public:
+  char*name;double*pval;
+  RVar(){name=NULL;pval=NULL;};
+  RVar(const RVar&);
+  RVar(const char*,double*);
+  ~RVar();
+  friend int operator==(const RVar&,const RVar&);
+};
+typedef RVar* PRVar;
+
+enum ROperator{ErrOp,Juxt,Num,Var,Add,Sub,Opp,Mult,Div,Pow,Sqrt,
+	       NthRoot,Abs,Sin,Cos,Tg,Ln,Exp,Acos,Asin,Atan,E10,Fun};
+
+typedef void ((*pfoncld)(double*&));
+
+class ROperation;
+typedef ROperation* PROperation;
+class RFunction;
+typedef RFunction* PRFunction;
+
+class ROperation{
+  pfoncld*pinstr;double**pvals;double*ppile;RFunction**pfuncpile;
+  mutable signed char containfuncflag;
+  void BuildCode();
+	void Destroy();
+ public:
+  ROperator op;
+  PROperation mmb1,mmb2;
+  double ValC;const RVar* pvar;double*pvarval;
+  RFunction* pfunc;
+  ROperation();
+  ROperation(const ROperation&);
+  ROperation(double);
+  ROperation(const RVar&);
+  ROperation(char*sp,int nvarp=0,PRVar*ppvarp=NULL,int nfuncp=0,PRFunction*ppfuncp=NULL);
+  ~ROperation();
+  double Val() const;
+  signed char ContainVar(const RVar&) const;
+  signed char ContainFunc(const RFunction&) const;
+  signed char ContainFuncNoRec(const RFunction&) const; // No recursive test on subfunctions
+  ROperation NthMember(int) const;int NMembers() const;
+  signed char HasError(const ROperation* =NULL) const;
+  ROperation& operator=(const ROperation&);
+  friend int operator==(const ROperation& ,const double);
+  friend int operator==(const ROperation& ,const ROperation&);
+  friend int operator!=(const ROperation& ,const ROperation&);
+  ROperation operator+() const;ROperation operator-() const;
+  friend ROperation operator,(const ROperation&,const ROperation&);
+  friend ROperation operator+(const ROperation&,const ROperation&);
+  friend ROperation operator-(const ROperation&,const ROperation&);
+  friend ROperation operator*(const ROperation&,const ROperation&);
+  friend ROperation operator/(const ROperation&,const ROperation&);
+  friend ROperation operator^(const ROperation&,const ROperation&);  // Caution: wrong associativity and precedence
+  friend ROperation sqrt(const ROperation&);
+  friend ROperation abs(const ROperation&);
+  friend ROperation sin(const ROperation&);
+  friend ROperation cos(const ROperation&);
+  friend ROperation tan(const ROperation&);
+  friend ROperation log(const ROperation&);
+  friend ROperation exp(const ROperation&);
+  friend ROperation acos(const ROperation&);
+  friend ROperation asin(const ROperation&);
+  friend ROperation atan(const ROperation&);
+  friend ROperation ApplyOperator(int,ROperation**,ROperation (*)(const ROperation&,const ROperation&));
+  ROperation Diff(const RVar&) const; //  Differentiate w.r.t a variable
+  char* Expr() const;
+  ROperation Substitute(const RVar&,const ROperation&) const;
+};
+
+class RFunction{
+  double*buf;
+public:
+  signed char type;
+  double ((*pfuncval)(double));
+  ROperation op;int nvars;RVar** ppvar;
+  char*name;
+  RFunction();
+  RFunction(double ((*)(double)));
+  RFunction(const ROperation& opp,RVar* pvarp);
+  RFunction(const ROperation& opp,int nvarsp,RVar**ppvarp);
+  RFunction(const RFunction&);
+  ~RFunction();
+  RFunction& operator=(const RFunction&);
+  void SetName(const char*s);
+  double Val(double) const;
+  double Val(double*) const;
+  friend int operator==(const RFunction&,const RFunction&);
+  ROperation operator()(const ROperation&);
+};
+
+
+char* MidStr(const char*s,int i1,int i2);
+char* CopyStr(const char*s);
+char* InsStr(const char*s,int n,char c);
+signed char EqStr(const char*s,const char*s2);
+signed char CompStr(const char*s,int n,const char*s2);
+char* DelStr(const char*s,int n);
+
+#endif
diff --git a/mathexpr/mathexpr/mathexpr_c.cpp b/mathexpr/mathexpr/mathexpr_c.cpp
new file mode 100644
index 0000000..7e7e906
--- /dev/null
+++ b/mathexpr/mathexpr/mathexpr_c.cpp
@@ -0,0 +1,1399 @@
+/*
+
+mathexpr_c.cpp version 2.0
+
+Copyright (c) 1997-2000 Yann OLLIVIER
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+*/
+
+#include"mathexpr_c.h"
+
+#ifndef __STD_COMPLEX_
+double_complex tan(double_complex z)
+{return sin(z)/cos(z);}
+double_complex acos(double_complex z)
+{return -double_complex(0,1)*log(z+double_complex(0,1)*sqrt(double_complex(1,0)-z*z));}
+double_complex asin(double_complex z)
+{return -double_complex(0,1)*log(double_complex(0,1)*z+sqrt(double_complex(1,0)-z*z));}
+double_complex atan(double_complex z)
+{return -double_complex(0,1)*log((double_complex(1,0)+double_complex(0,1)*z)/sqrt(double_complex(1,0)+z*z));}
+#endif
+
+char* MidStr(const char*s,int i1,int i2)
+{
+ if(i1<0||i2>=(int)strlen(s)||i1>i2){
+   char* cp = new char[1];
+   cp[0] = '\0';
+   return cp;
+ }
+ char*s1=new char[i2-i1+2];
+ int i;
+ for(i=i1;i<=i2;i++)s1[i-i1]=s[i];
+ s1[i2-i1+1]=0;return s1;
+}
+
+char* CopyStr(const char*s)
+{char*s1=new char[strlen(s)+1];char*s12=s1;const char*s2=s;
+while((*s12++=*s2++));return s1;}
+
+void InsStr(char*&s,int n,char c)// Warning : deletes the old string
+{if(n<0||n>(int)strlen(s))return;
+char*s1=new char[strlen(s)+2];
+int i;
+for(i=0;i<n;i++)s1[i]=s[i];
+s1[n]=c;for(i=n+1;s[i-1];i++)s1[i]=s[i-1];
+s1[i]=0;
+delete[]s;s=s1;
+}
+
+signed char EqStr(const char*s,const char*s2)
+{if(strlen(s)!=strlen(s2))return 0;
+int i;
+for(i=0;i<s[i];i++)if(s[i]!=s2[i])return 0;
+return 1;
+}
+
+signed char CompStr(const char*s,int n,const char*s2)
+{if(n<0||n>=(int)strlen(s)||n+(int)strlen(s2)>(int)strlen(s))return 0;
+int i;
+for(i=0;s2[i];i++)if(s[i+n]!=s2[i])return 0;
+return 1;
+}
+
+void DelStr(char*&s,int n)//Deletes the old string
+{char*s1=new char[strlen(s)];
+int i;
+for(i=0;i<n;i++)s1[i]=s[i];
+for(i=n;s[i+1];i++)s1[i]=s[i+1];
+s1[i]=0;
+delete[]s;s=s1;
+}
+
+CVar::CVar(const CVar & rvarp)
+{if(this==&rvarp)return;pval=rvarp.pval;name=CopyStr(rvarp.name);
+}
+
+CVar::CVar(const char*namep,double_complex*pvalp)
+{pval=(double_complex*)pvalp;name=CopyStr(namep);}
+
+CVar::~CVar()
+{if(name!=NULL)delete[] name;}
+
+CFunction::CFunction()
+{
+  type=-1;name=new char[1];name[0]=0;
+  nvars=0;ppvar=NULL;pfuncval=NULL;op=ErrVal;buf=NULL;
+}
+
+CFunction::CFunction(double_complex ((*pfuncvalp)(double_complex)))
+{
+  type=0;pfuncval=pfuncvalp;name=new char[1];name[0]=0;
+  nvars=1;ppvar=NULL;op=ErrVal;buf=NULL;
+}
+
+CFunction::CFunction(const CFunction& rfunc)
+{
+  if(this==&rfunc)return;
+  type=rfunc.type;op=rfunc.op;
+  pfuncval=rfunc.pfuncval;
+  name=CopyStr(rfunc.name);
+  nvars=rfunc.nvars;
+  if(rfunc.ppvar!=NULL&&nvars){
+    ppvar=new PCVar[nvars];
+    int i;for(i=0;i<nvars;i++)ppvar[i]=rfunc.ppvar[i];
+    buf=new double_complex[nvars];
+  }else {ppvar=NULL;buf=NULL;}
+}
+
+CFunction::CFunction(const COperation& opp,CVar* pvarp):op(opp)
+{
+  type=1;name=new char[1];name[0]=0;
+  nvars=1;ppvar=new PCVar[1];ppvar[0]=pvarp;buf=new double_complex[1];
+}
+
+CFunction::CFunction(const COperation& opp, int nvarsp,CVar**ppvarp):op(opp)
+{
+  type=1;name=new char[1];name[0]=0;
+  nvars=nvarsp;
+  if(nvars){
+    ppvar=new PCVar[nvars];
+    int i;for(i=0;i<nvars;i++)ppvar[i]=ppvarp[i];
+    buf=new double_complex[nvars];
+  }else {ppvar=NULL;buf=NULL;}
+}
+
+CFunction::~CFunction()
+{
+  if(name!=NULL)delete[]name;
+  if(ppvar!=NULL)delete[]ppvar;
+  if(buf!=NULL)delete[]buf;
+}
+
+CFunction& CFunction::operator=(const CFunction& rfunc)
+{
+  if(this==&rfunc)return *this;
+  type=rfunc.type;op=rfunc.op;
+  pfuncval=rfunc.pfuncval;
+  delete[]name;
+  name=CopyStr(rfunc.name);
+  if(ppvar!=NULL)delete[]ppvar;
+  ppvar=NULL;
+  if(buf!=NULL)delete[]buf;
+  buf=NULL;
+  nvars=rfunc.nvars;
+  if(type==1&&nvars){
+    ppvar=new PCVar[nvars];buf=new double_complex[nvars];
+    int i;for(i=0;i<nvars;i++)ppvar[i]=rfunc.ppvar[i];
+  }
+  return *this;
+}
+
+void CFunction::SetName(const char*s)
+{if(name!=NULL)delete[]name;name=CopyStr(s);}
+
+double_complex CFunction::Val(double_complex x) const
+{
+  if(type==-1||nvars>=2)return ErrVal;
+  if(type==0)return (*pfuncval)(x);
+  double_complex xb=*(*ppvar)->pval,y;
+  *(*ppvar)->pval=x;  // Warning : could cause trouble if this value is used in a parallel process
+  y=op.Val();
+  *(*ppvar)->pval=xb;
+  return y;
+}
+
+double_complex CFunction::Val(double_complex*pv) const
+{
+  if(type==-1)return ErrVal;
+  if(type==0)return (*pfuncval)(*pv);
+  double_complex y;
+  int i;
+  for(i=0;i<nvars;i++){
+    buf[i]=*ppvar[i]->pval;
+    // Warning : could cause trouble if this value is used in a parallel process
+    *ppvar[i]->pval=pv[i];
+  }
+  y=op.Val();
+  for(i=0;i<nvars;i++)*ppvar[i]->pval=buf[i];
+  return y;
+}
+
+COperation::COperation()
+{op=ErrOp;mmb1=NULL;mmb2=NULL;ValC=ErrVal;pvar=NULL;pvarval=NULL;pfunc=NULL;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;BuildCode();}
+
+COperation::~COperation()
+{
+  Destroy();
+}
+
+COperation::COperation(const COperation&COp)
+{
+  op=COp.op;pvar=COp.pvar;pvarval=COp.pvarval;ValC=COp.ValC;pfunc=COp.pfunc;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  if(COp.mmb1!=NULL)mmb1=new COperation(*(COp.mmb1));else mmb1=NULL;
+  if(COp.mmb2!=NULL)mmb2=new COperation(*(COp.mmb2));else mmb2=NULL;
+  BuildCode();
+}
+
+COperation::COperation(double_complex x)
+{
+  if(x==ErrVal){op=ErrOp;mmb1=NULL;mmb2=NULL;ValC=ErrVal;}
+  else {op=Num;mmb1=NULL;mmb2=NULL;ValC=x;}
+  pvar=NULL;pvarval=NULL;pfunc=NULL;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  BuildCode();
+}
+
+COperation::COperation(double x)
+{
+  if(x==ErrVal){op=ErrOp;mmb1=NULL;mmb2=NULL;ValC=ErrVal;}
+  else if(x>=0){op=Num;mmb1=NULL;mmb2=NULL;ValC=x;}
+  else{op=Opp;mmb1=NULL;mmb2=new COperation(-x);ValC=ErrVal;}
+  pvar=NULL;pvarval=NULL;pfunc=NULL;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  BuildCode();
+}
+
+COperation::COperation(const CVar& varp)
+{op=Var;mmb1=NULL;mmb2=NULL;ValC=ErrVal;pvar=&varp;pvarval=varp.pval;containfuncflag=0;pfunc=NULL;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;BuildCode();}
+
+COperation& COperation::operator=(const COperation& COp)
+{
+  if(this==&COp)return *this;
+  Destroy();
+  op=COp.op;pvar=COp.pvar;pvarval=COp.pvarval;ValC=COp.ValC;pfunc=COp.pfunc;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  if(COp.mmb1!=NULL)mmb1=new COperation(*(COp.mmb1));else mmb1=NULL;
+  if(COp.mmb2!=NULL)mmb2=new COperation(*(COp.mmb2));else mmb2=NULL;
+  BuildCode();
+  return *this;
+}
+
+int operator==(const COperation& op,const double v)
+{return(op.op==Num&&op.ValC==v);}
+
+int operator==(const COperation& op,const double_complex v)
+{return(op.op==Num&&op.ValC==v);}
+
+int operator==(const COperation& op1,const COperation& op2)
+{if(op1.op!=op2.op)return 0;
+if(op1.op==Var)return(*(op1.pvar)==*(op2.pvar));
+if(op1.op==Fun)return(op1.pfunc==op2.pfunc); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
+if(op1.op==Num)return(op1.ValC==op2.ValC);
+if(op1.mmb1==NULL&&op2.mmb1!=NULL)return 0;
+if(op1.mmb2==NULL&&op2.mmb2!=NULL)return 0;
+if(op2.mmb1==NULL&&op1.mmb1!=NULL)return 0;
+if(op2.mmb2==NULL&&op1.mmb2!=NULL)return 0;
+return(((op1.mmb1==NULL&&op2.mmb1==NULL)||(*(op1.mmb1)==*(op2.mmb1)))&&
+	((op1.mmb2==NULL&&op2.mmb2==NULL)||(*(op1.mmb2)==*(op2.mmb2))));
+}
+
+int operator!=(const COperation& op1,const COperation& op2)
+{
+  if(op1.op!=op2.op)return 1;
+  if(op1.op==Var)return(op1.pvar!=op2.pvar);
+  if(op1.op==Fun)return(!(op1.pfunc==op2.pfunc)); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
+  if(op1.op==Num)return(op1.ValC!=op2.ValC);
+  if(op1.mmb1==NULL&&op2.mmb1!=NULL)return 1;
+  if(op1.mmb2==NULL&&op2.mmb2!=NULL)return 1;
+  if(op2.mmb1==NULL&&op1.mmb1!=NULL)return 1;
+  if(op2.mmb2==NULL&&op1.mmb2!=NULL)return 1;
+  return(((op1.mmb1!=NULL||op2.mmb1!=NULL)&&(*(op1.mmb1)!=*(op2.mmb1)))||
+	 ((op1.mmb2!=NULL||op2.mmb2!=NULL)&&(*(op1.mmb2)!=*(op2.mmb2))));
+}
+
+COperation COperation::operator+() const
+{return *this;}
+
+COperation COperation::operator-() const
+{if(op==Num)return -ValC;
+COperation resultat;
+if(op==Opp)resultat=*mmb2;else{resultat.op=Opp;resultat.mmb2=new COperation(*this);};
+return resultat;
+}
+
+COperation operator,(const COperation& op1,const COperation& op2)
+{COperation resultat;
+resultat.op=Juxt;resultat.mmb1=new COperation(op1);
+resultat.mmb2=new COperation(op2);
+return resultat;
+}
+
+COperation operator+(const COperation& op1,const COperation& op2)
+{
+if(op1.op==Num&&op2.op==Num)return op1.ValC+op2.ValC;
+if(op1==0.)return op2;if(op2==0.)return op1;
+if(op1.op==Opp)return op2-*(op1.mmb2);if(op2.op==Opp)return op1-*(op2.mmb2);
+COperation resultat;
+resultat.op=Add;resultat.mmb1=new COperation(op1);
+resultat.mmb2=new COperation(op2);
+return resultat;
+}
+
+COperation operator-(const COperation& op1,const COperation& op2)
+{
+if(op1.op==Num&&op2.op==Num)return op1.ValC-op2.ValC;
+if(op1==0.)return -op2;if(op2==0.)return op1;
+if(op1.op==Opp)return -(op2+*(op1.mmb2));if(op2.op==Opp)return op1+*(op2.mmb2);
+COperation resultat;
+resultat.op=Sub;resultat.mmb1=new COperation(op1);
+resultat.mmb2=new COperation(op2);
+return resultat;
+}
+
+COperation operator*(const COperation& op1,const COperation& op2)
+{
+if(op1.op==Num&&op2.op==Num)return op1.ValC*op2.ValC;
+if(op1==0.||op2==0.)return 0.;
+if(op1==1.)return op2;if(op2==1.)return op1;
+if(op1.op==Opp)return -(*(op1.mmb2)*op2);if(op2.op==Opp)return -(op1**(op2.mmb2));
+COperation resultat;
+resultat.op=Mult;resultat.mmb1=new COperation(op1);
+resultat.mmb2=new COperation(op2);
+return resultat;
+}
+
+COperation operator/(const COperation& op1,const COperation& op2)
+{if(op1.op==Num&&op2.op==Num)return (op2.ValC!=double_complex(0,0)?op1.ValC/op2.ValC:ErrVal);
+if(op1==0.0)return 0.;if(op2==1.)return op1;if(op2==0.)return ErrVal;
+if(op1.op==Opp)return -(*(op1.mmb2)/op2);if(op2.op==Opp)return -(op1/(*(op2.mmb2)));
+COperation resultat;
+resultat.op=Div;resultat.mmb1=new COperation(op1);
+resultat.mmb2=new COperation(op2);
+return resultat;
+}
+
+COperation operator^(const COperation& op1,const COperation& op2)
+{if(op1==0.)return 0.;
+if(op2==0.)return 1.;
+if(op2==1.)return op1;
+COperation resultat;
+resultat.op=Pow;resultat.mmb1=new COperation(op1);
+resultat.mmb2=new COperation(op2);
+return resultat;
+}               
+                
+COperation real(const COperation& op)
+{
+  if(op.op==Num)return real(op.ValC);
+  COperation rop;rop.op=Real;rop.mmb2=new COperation(op);return rop;
+}
+COperation imag(const COperation& op)
+{
+  if(op.op==Num)return imag(op.ValC);
+  COperation rop;rop.op=Imag;rop.mmb2=new COperation(op);return rop;
+}
+COperation conj(const COperation& op)
+{
+  if(op.op==Num)return conj(op.ValC);
+  COperation rop;rop.op=Conj;rop.mmb2=new COperation(op);return rop;
+}
+COperation arg(const COperation& op)
+{COperation rop;rop.op=Arg;rop.mmb2=new COperation(op);return rop;}
+COperation sqrt(const COperation& op)
+{COperation rop;rop.op=Sqrt;rop.mmb2=new COperation(op);return rop;}
+COperation abs(const COperation& op)
+{COperation rop;rop.op=Abs;rop.mmb2=new COperation(op);return rop;}
+COperation sin(const COperation& op)
+{COperation rop;rop.op=Sin;rop.mmb2=new COperation(op);return rop;}
+COperation cos(const COperation& op)
+{COperation rop;rop.op=Cos;rop.mmb2=new COperation(op);return rop;}
+COperation tan(const COperation& op)
+{COperation rop;rop.op=Tg;rop.mmb2=new COperation(op);return rop;}
+COperation log(const COperation& op)
+{COperation rop;rop.op=Ln;rop.mmb2=new COperation(op);return rop;}
+COperation exp(const COperation& op)
+{COperation rop;rop.op=Exp;rop.mmb2=new COperation(op);return rop;}
+COperation acos(const COperation& op)
+{COperation rop;rop.op=Acos;rop.mmb2=new COperation(op);return rop;}
+COperation asin(const COperation& op)
+{COperation rop;rop.op=Asin;rop.mmb2=new COperation(op);return rop;}
+COperation atan(const COperation& op)
+{COperation rop;rop.op=Atan;rop.mmb2=new COperation(op);return rop;}
+
+COperation ApplyOperator(int n,COperation**pops,COperation (*func)(const COperation&,const COperation&))
+{
+  if(n<=0)return ErrVal;
+  if(n==1)return *pops[0];
+  if(n==2)return (*func)(*pops[0],*pops[1]);
+  return (*func)(*pops[0],ApplyOperator(n-1,pops+1,func));
+}
+
+COperation CFunction::operator()(const COperation& op)
+{
+  /* Code to use to replace explcitly instead of using a pointer to
+     if(nvars!=op.NMembers()||type==-1||type==0)return ErrVal;
+     COperation op2=*pop;int i;
+     CVar**ppvar2=new PCVar[nvars];char s[11]="";
+     for(i=0;i<nvars;i++){
+     sprintf(s,";var%i;",i);
+     ppvar2[i]=new CVar(s,NULL);
+     op2=op2.Substitute(*ppvar[i],(COperation)*ppvar2[i]);
+     }
+     for(i=0;i<nvars;i++){
+     op2=op2.Substitute(*ppvar2[i],op.NthMember(i+1));
+     delete ppvar2[i];
+     }
+     delete[]ppvar2;
+     return op2;
+  */
+  COperation op2;op2.op=Fun;op2.pfunc=this;op2.mmb2=new COperation(op);
+  return op2;
+}
+
+//Auxiliary string functions
+
+void SupprSpaces(char*&s)//Deletes the old string
+{
+ int i;
+ for(i=0;s[i];i++)if(s[i]==' '||s[i]=='\t'||s[i]=='\n')DelStr(s,i--);
+}
+
+signed char IsNumeric(char c)
+{if(c!='0'&&c!='1'&&c!='2'&&c!='3'&&c!='4'
+	&&c!='5'&&c!='6'&&c!='7'&&c!='8'&&c!='9'&&c!='.')return 0;
+return 1;
+}
+
+signed char IsTNumeric(char *s)
+{int i;for(i=0;i<(int)strlen(s);i++)if(!IsNumeric(s[i]))return 0;return 1;
+}
+
+int SearchCorOpenbracket(char*s,int n)  //Searchs the corresponding bracket of an opening bracket
+{if(n>=(int)strlen(s)-1)return -1;
+int i,c=1;
+for(i=n+1;s[i];i++){
+	if(s[i]=='(')c++;else if(s[i]==')')c--;
+	if(!c)return i;
+  };
+return -1;
+}
+
+int SearchCorClosebracket(char*s,int n)  //Searchs the corresponding bracket of a closing bracket
+{if(n<1)return -1;
+int i,c=1;
+for(i=n-1;i>=0;i--){
+	if(s[i]==')')c++;else if(s[i]=='(')c--;
+	if(!c)return i;
+  };
+return -1;
+}
+
+int SearchOperator(char*s,COperator op)
+{
+  char opc;
+  switch(op){
+  case ErrOp:case Num:case Var:return -1;
+  case Juxt:opc=',';break;
+  case Add:opc='+';break;
+  case Sub:opc='-';break;
+  case Mult:opc='*';break;
+  case Div:opc='/';break;
+  case Pow:opc='^';break;
+  case NthRoot:opc='#';break;
+  case E10:opc='E';break;        
+  default:return -1;
+  };
+  int i;
+  for(i=(int)strlen(s)-1;i>=0;i--){
+    if(s[i]==opc&&(op!=Sub||i&&s[i-1]==')'))return i;
+    if(s[i]==')'){i=SearchCorClosebracket(s,i);if(i==-1)return -1;};
+  };
+  return -1;
+}
+
+void SimplifyStr(char*&s) //Warning : deletes the old string
+{if(!strlen(s))return;
+char*s1=s,*s2=s+strlen(s);signed char ind=0;
+if(s1[0]=='('&&SearchCorOpenbracket(s1,0)==s2-s1-1){
+  s1++;s2--;ind=1;}
+if(s1==s2)
+{
+	delete[]s;
+	s=new char[1]; // ISO C++ forbids initialization in array new
+	s[0]=0;
+	return;
+}
+if(s1[0]==' '){ind=1;while(s1[0]==' '&&s1<s2)s1++;}
+if(s1==s2)
+{
+	delete[]s;
+	s=new char[1]; // ISO C++ forbids initialization in array new
+	s[0]=0;
+	return;
+}
+if(*(s2-1)==' '){ind=1;while(s2>s1&&*(s2-1)==' ')s2--;}
+*s2=0;
+s1=CopyStr(s1);delete[]s;s=s1;
+if(ind)SimplifyStr(s);
+}
+
+int max(int a, int b){return (a>b?a:b);}
+
+int IsVar(const char*s,int n,int nvar,PCVar*ppvar)
+{
+  if(n<0||n>(int)strlen(s))return 0;
+  int i;int l=0;
+  for(i=0;i<nvar;i++)if(CompStr(s,n,(*(ppvar+i))->name))l=max(l,strlen((*(ppvar+i))->name));
+  return l;
+}
+
+int IsFunction(const char*s,int n)
+{
+  if(CompStr(s,n,"sin")||CompStr(s,n,"cos")||CompStr(s,n,"exp")
+     ||CompStr(s,n,"tan")||CompStr(s,n,"log")||CompStr(s,n,"atg")
+     ||CompStr(s,n,"arg")||CompStr(s,n,"abs"))return 3;
+  if(CompStr(s,n,"tg")||CompStr(s,n,"ln"))return 2;
+  if(CompStr(s,n,"sqrt")||CompStr(s,n,"asin")||CompStr(s,n,"atan")
+     ||CompStr(s,n,"real")||CompStr(s,n,"conj")
+     ||CompStr(s,n,"imag")||CompStr(s,n,"acos"))return 4;
+  if(CompStr(s,n,"arcsin")||CompStr(s,n,"arccos")||CompStr(s,n,"arctan"))return 6;
+  if(CompStr(s,n,"arctg"))return 5;
+  return 0;
+}
+
+int IsFunction(const char*s,int n,int nfunc,PCFunction*ppfunc)
+  //Not recognized if a user-defined function is eg "sine" ie begins like
+  //a standard function
+  //IF PATCHED TO DO OTHERWISE, SHOULD BE PATCHED TOGETHER WITH THE
+  //PARSER BELOW which treats standard functions before user-defined ones
+{
+  int l=IsFunction(s,n);
+  if(l)return l;
+  int i;l=0;
+  for(i=0;i<nfunc;i++)if(CompStr(s,n,ppfunc[i]->name))l=max(l,strlen(ppfunc[i]->name));
+  return l;
+}
+
+signed char IsFunction(COperator op)
+{return (op==Exp||op==Abs||op==Sin||op==Cos||op==Tg||op==Ln
+	||op==Atan||op==Asin||op==Acos||op==Atan||op==Sqrt||op==Opp
+	||op==Real||op==Imag||op==Conj||op==Arg);
+}
+
+void IsolateVars(char*&s,int nvar,PCVar*ppvar,int nfunc,PCFunction*ppfunc)//Deletes the old string
+{
+  int i,j;
+  i=0;
+  for(i=0;s[i];i++){
+    if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)return;continue;};
+    if(((j=IsVar(s,i,nvar,ppvar))>IsFunction(s,i,nfunc,ppfunc))||((CompStr(s,i,"pi")||CompStr(s,i,"PI")||CompStr(s,i,"Pi"))&&(j=2)
+       ||(CompStr(s,i,"i")&&(j=1)))){
+      InsStr(s,i,'(');InsStr(s,i+j+1,')');i+=j+1;continue;};
+    if(IsFunction(s,i,nfunc,ppfunc)){i+=IsFunction(s,i,nfunc,ppfunc)-1;if(!s[i])return;continue;};
+  };
+}
+
+void IsolateNumbers(char*&s,int nvar,CVar**ppvar,int nfunc,CFunction**ppfunc)//Deletes the old string
+{
+  int i,i2,ind=0,t1,t2;
+  for(i=0;s[i];i++){
+    if(ind&&!IsNumeric(s[i])){ind=0;InsStr(s,i2,'(');i++;InsStr(s,i,')');continue;};
+		t1=IsVar(s,i,nvar,ppvar);t2=IsFunction(s,i,nfunc,ppfunc);
+		if(t1||t2){i+=max(t1,t2)-1;continue;}
+    if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)return;continue;};
+    if(!ind&&IsNumeric(s[i])){i2=i;ind=1;};
+  };
+if(ind)InsStr(s,i2,'(');i++;InsStr(s,i,')');
+}
+
+COperation::COperation(char*sp,int nvar,PCVar*ppvarp,int nfuncp,PCFunction*ppfuncp)
+{
+  ValC=ErrVal;mmb1=NULL;mmb2=NULL;pvar=NULL;op=ErrOp;pvarval=NULL;containfuncflag=0;pfunc=NULL;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+  int i,j,k,l;signed char flag=1;
+  char*s=CopyStr(sp),*s1=NULL,*s2=NULL;
+  SimplifyStr(s);if(!s[0]||!strcmp(s,"Error")){goto fin;}
+  while(s[0]==':'||s[0]==';'){
+		s1=CopyStr(s+1);delete[]s;s=s1;s1=NULL;
+    SimplifyStr(s);if(!s[0]||!strcmp(s,"Error")){goto fin;}
+  }
+  if(IsTNumeric(s)){op=Num;ValC=atof(s);mmb1=NULL;mmb2=NULL;goto fin;};
+  if(EqStr(s,"pi")||EqStr(s,"PI")||EqStr(s,"Pi"))
+    {op=Num;ValC=3.141592653589793238462643383279L;mmb1=NULL;mmb2=NULL;goto fin;};
+  if(EqStr(s,"i")||EqStr(s,"I"))
+    {op=Num;ValC=double_complex(0,1);mmb1=NULL;mmb2=NULL;goto fin;};
+  if(IsFunction(s,0,nfuncp,ppfuncp)<IsVar(s,0,nvar,ppvarp))
+    for(i=0;i<nvar;i++)if(EqStr(s,(*(ppvarp+i))->name))
+      {pvar=ppvarp[i];pvarval=pvar->pval;op=Var;mmb1=NULL;mmb2=NULL;goto fin;};
+  for(k=0;s[k];k++){
+    if(s[k]=='('){k=SearchCorOpenbracket(s,k);if(k==-1)break;continue;};
+    if((l=IsFunction(s,k,nfuncp,ppfuncp))&&l>=IsVar(s,k,nvar,ppvarp)){
+      i=k+l;while(s[i]==' ')i++;
+      if(s[i]=='('){
+				j=SearchCorOpenbracket(s,i);
+				if(j!=-1){InsStr(s,i,';');k=j+1;}
+			}else if(s[i]!=':'&&s[i]!=';'){InsStr(s,i,':');k=i;}
+    }
+  }
+  IsolateNumbers(s,nvar,ppvarp,nfuncp,ppfuncp);
+	if(nvar)IsolateVars(s,nvar,ppvarp,nfuncp,ppfuncp);
+  SupprSpaces(s);
+  i=SearchOperator(s,Juxt);
+  if(i!=-1){char*s1="",*s2="";s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Juxt;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Add);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Add;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Sub);
+  if(i!=-1){
+    s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+    op=Sub;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+    mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  if(s[0]=='-'){s2=MidStr(s,1,strlen(s)-1);
+  op=Opp;mmb1=NULL;
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  for(i=0;s[i];i++){
+		if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)break;continue;};
+		if(IsFunction(s,i,nfuncp,ppfuncp)){
+			k=i+IsFunction(s,i,nfuncp,ppfuncp);while(s[k]==' ')k++;
+			if(s[k]==';'){
+				//	s=DelStr(s,k);
+				j=k;while(s[j]!='(')j++;
+				j=SearchCorOpenbracket(s,j);
+				if(j!=-1){InsStr(s,j,')');InsStr(s,i,'(');i=j+2;}
+			}else if(s[k]==':'){
+				//	s=DelStr(s,k);
+				for(j=k;s[j];j++)
+					if(s[j]=='('){j=SearchCorOpenbracket(s,j);break;}
+				if(j==-1)break;
+				for(j++;s[j];j++){
+					if(s[j]=='('){j=SearchCorOpenbracket(s,j);if(j==-1){flag=0;break;};continue;};
+					if(IsFunction(s,j,nfuncp,ppfuncp))break;
+				}
+				if(flag==0){flag=1;break;}
+				while(j>i&&s[j-1]!=')')j--;if(j<=i+1)break;
+				InsStr(s,i,'(');InsStr(s,j+1,')');
+				i=j+1;
+			}
+		}
+  }
+  for(i=0;s[i]&&s[i+1];i++)if(s[i]==')'&&s[i+1]=='(')
+    InsStr(s,++i,'*');
+  if(s[0]=='('&&SearchCorOpenbracket(s,0)==(int)strlen(s)-1){
+    if(CompStr(s,1,"exp")){op=Exp;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"real")){op=Real;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"imag")){op=Imag;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"conj")){op=Conj;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"arg")){op=Arg;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"abs")){op=Abs;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"sin")){op=Sin;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"cos")){op=Cos;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"tan")){op=Tg;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"log")){op=Ln;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"atg")){op=Atan;s2=MidStr(s,4,strlen(s)-2);}
+    else if(CompStr(s,1,"tg")){op=Tg;s2=MidStr(s,3,strlen(s)-2);}
+    else if(CompStr(s,1,"ln")){op=Ln;s2=MidStr(s,3,strlen(s)-2);}
+    else if(CompStr(s,1,"asin")){op=Asin;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"acos")){op=Acos;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"atan")){op=Atan;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"sqrt")){op=Sqrt;s2=MidStr(s,5,strlen(s)-2);}
+    else if(CompStr(s,1,"arcsin")){op=Asin;s2=MidStr(s,7,strlen(s)-2);}
+    else if(CompStr(s,1,"arccos")){op=Acos;s2=MidStr(s,7,strlen(s)-2);}
+    else if(CompStr(s,1,"arctan")){op=Atan;s2=MidStr(s,7,strlen(s)-2);}
+    else if(CompStr(s,1,"arctg")){op=Atan;s2=MidStr(s,6,strlen(s)-2);}
+    else {
+      for(i=-1,k=0,j=0;j<nfuncp;j++)if(CompStr(s,1,ppfuncp[j]->name)&&k<(int)strlen(ppfuncp[j]->name)){k=strlen(ppfuncp[j]->name);i=j;}
+      if(i>-1){
+				op=Fun;s2=MidStr(s,strlen(ppfuncp[i]->name)+1,strlen(s)-2);
+				pfunc=ppfuncp[i];
+      }
+    }
+    mmb1=NULL;mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);
+    if(op==Fun)if(mmb2->NMembers()!=pfunc->nvars){op=ErrOp;mmb1=NULL;mmb2=NULL;goto fin;}
+    goto fin;
+  };
+  i=SearchOperator(s,Mult);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Mult;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Div);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Div;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,Pow);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=Pow;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,NthRoot);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  if(i==0||s[i-1]!=')')
+    {op=Sqrt;mmb1=NULL;}else
+      {op=NthRoot;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);};
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  i=SearchOperator(s,E10);
+  if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+  op=E10;mmb1=new COperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+  mmb2=new COperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
+  };
+  op=ErrOp;mmb1=NULL;mmb2=NULL;
+fin:
+	BuildCode();
+	delete[]s;if(s1!=NULL)delete[] s1;if(s2!=NULL)delete[]s2;
+}
+
+void COperation::Destroy()
+{
+  if(mmb1!=NULL&&mmb2!=NULL&&mmb1!=mmb2){delete mmb1;delete mmb2;mmb1=NULL;mmb2=NULL;}
+  else if(mmb1!=NULL){delete mmb1;mmb1=NULL;}else if(mmb2!=NULL){delete mmb2;mmb2=NULL;}
+  if(pinstr!=NULL){delete[]pinstr;pinstr=NULL;}
+	if(pvals!=NULL){
+		if(op==ErrOp||op==Num)delete pvals[0];
+		delete[]pvals;pvals=NULL;
+	}
+  if(ppile!=NULL){delete[]ppile;ppile=NULL;}
+  if(pfuncpile!=NULL){delete[]pfuncpile;pfuncpile=NULL;}
+}
+
+int operator==(const CVar& var1,const CVar& var2)
+{return(var1.pval==var2.pval&&EqStr(var1.name,var2.name));
+}
+
+int operator==(const CFunction& f1,const CFunction& f2)
+{
+  if(f1.type!=f2.type)return 0;
+  if(f1.type==-1)return 1; // Nonfunction==nonfunction
+  if(f1.type==0)return (f1.pfuncval==f2.pfuncval&&EqStr(f1.name,f2.name));
+  if(f1.op!=f2.op)return 0;
+  if(!EqStr(f1.name,f2.name))return 0;
+  if(f1.nvars!=f2.nvars)return 0;
+  int i;
+  for(i=0;i<f1.nvars;i++)if(!(*f1.ppvar[i]==*f2.ppvar[i]))return 0;
+  return 1;
+}
+
+/*
+double_complex COperation::Val() const // Won't work if multi-variable functions are included
+{
+  double_complex v1=ErrVal,v2=ErrVal;
+  if(mmb1!=NULL){v1=mmb1->Val();if(norm(v1)<sqrtminfloat)v1=0;else if(v1==double_complex(ErrVal)||norm(v1)>sqrtmaxfloat)return ErrVal;};
+  if(mmb2!=NULL){v2=mmb2->Val();if(norm(v2)<sqrtminfloat)v2=0;else if(v2==double_complex(ErrVal)||norm(v2)>sqrtmaxfloat)return ErrVal;};
+  switch(op){
+  case Num:return ValC;
+  case Var:return *pvarval;
+  case Add:return v1+v2;
+  case Sub:return v1-v2;
+  case Opp:return -v2;
+  case Mult:return v1*v2;
+  case Div:if(v2!=double_complex(0))return v1/v2;else return ErrVal;
+  case Pow:if(v1==double_complex(0))return 0;else
+    if(norm(v2)*log(norm(v1))<DBL_MAX_EXP)return pow(v1,v2);else return ErrVal;
+  case Sqrt:return sqrt(v2);
+  case NthRoot:if(v1==double_complex(0)||norm(v2)*log(fabs(norm(v1)))<DBL_MIN_EXP)return ErrVal;
+  else return pow(v2,1/v1);
+  case E10:if(norm(v2)<300)return v1*pow(double_complex(10.),v2);else return ErrVal;
+  case Ln:if(v2!=double_complex(0))return log(v2);else return ErrVal;
+  case Exp:if(real(v2)<DBL_MAX_EXP)return exp(v2);else return ErrVal;
+  case Sin:if(norm(v2)<DBL_MAX_EXP)
+	    return sin(v2);else return ErrVal;
+  case Cos:if(norm(v2)<DBL_MAX_EXP)return cos(v2);else return ErrVal;
+  case Tg:if(norm(v2)<DBL_MAX_EXP)
+	   return tan(v2);else return ErrVal;
+  case Atan:if(v2==double_complex(0,1)||v2==double_complex(0,-1))return ErrVal;
+  else return atan(v2);
+  case Asin:return asin(v2);
+  case Acos:return acos(v2);
+  case Abs:return norm(v2);
+  case Real:return real(v2);
+  case Imag:return imag(v2);
+  case Conj:return conj(v2);
+  case Arg:if(v2!=double_complex(0))return arg(v2);else return ErrVal; 
+  case Fun:return pfunc->Val(v2);
+  default:return ErrVal;
+  };
+}
+*/
+
+signed char COperation::ContainVar(const CVar& varp) const
+{if(op==Var){if(EqStr(pvar->name,varp.name)&&pvar->pval==varp.pval)
+	return 1;else return 0;};
+if(mmb1!=NULL&&mmb1->ContainVar(varp))return 1;
+if(mmb2!=NULL&&mmb2->ContainVar(varp))return 1;
+return 0;
+}
+
+signed char COperation::ContainFuncNoRec(const CFunction& func) const // No recursive test on subfunctions
+{if(op==Fun){if(*pfunc==func)
+  return 1;else return 0;}
+ if(mmb1!=NULL&&mmb1->ContainFuncNoRec(func))return 1;
+ if(mmb2!=NULL&&mmb2->ContainFuncNoRec(func))return 1;
+ return 0;
+}
+
+signed char COperation::ContainFunc(const CFunction& func) const // Recursive test on subfunctions
+{
+  if(containfuncflag)return 0;
+  if(op==Fun&&*pfunc==func)return 1;
+  containfuncflag=1;
+  if(op==Fun)if(pfunc->op.ContainFunc(func)){containfuncflag=0;return 1;}
+  if(mmb1!=NULL&&mmb1->ContainFunc(func)){containfuncflag=0;return 1;}
+  if(mmb2!=NULL&&mmb2->ContainFunc(func)){containfuncflag=0;return 1;}
+  containfuncflag=0;return 0;
+}
+
+signed char COperation::HasError(const COperation*pop) const
+{
+  if(op==ErrOp)return 1;
+  if(op==Fun&&pfunc->type==1&&pfunc->op==*(pop==NULL?this:pop))return 1;
+  if(op==Fun&&pfunc->type==1&&pfunc->op.HasError((pop==NULL?this:pop)))return 1;
+  if(mmb1!=NULL&&mmb1->HasError((pop==NULL?this:pop)))return 1;
+  if(mmb2!=NULL&&mmb2->HasError((pop==NULL?this:pop)))return 1;
+  if(op==Fun&&pfunc->type==-1)return 1;
+  return 0;
+}
+
+int COperation::NMembers() const //Number of members for an operation like a,b,c...
+{
+  if(op==Fun)return(pfunc->type==1?pfunc->op.NMembers():pfunc->type==0?1:0);
+  if(op!=Juxt)return 1;else if(mmb2==NULL)return 0;else return 1+mmb2->NMembers();
+}
+
+COperation COperation::NthMember(int n) const
+{
+  PCFunction prf;
+  if(op==Fun&&pfunc->type==1&&pfunc->op.NMembers()>1){
+    prf=new CFunction(pfunc->op.NthMember(n),pfunc->nvars,pfunc->ppvar);
+    char*s=new char[strlen(pfunc->name)+10];
+    sprintf(s,"(%s_%i)",pfunc->name,n);prf->SetName(s);delete[]s;
+    return(*prf)(*mmb2);
+  }
+  if(n==1){
+    if(op!=Juxt)return *this; else if(mmb1!=NULL)return *mmb1;else return ErrVal;
+  };
+  if(op!=Juxt)return ErrVal;
+  if(n>1&&mmb2!=NULL)return mmb2->NthMember(n-1);
+  return ErrVal;
+}
+
+COperation COperation::Substitute(const CVar& var,const COperation& rop) const // Replaces variable var with expression rop
+{
+  if(!ContainVar(var))return *this;
+  if(op==Var)return rop;
+  COperation r;
+  r.op=op;r.pvar=pvar;r.pvarval=pvarval;r.ValC=ValC;r.pfunc=pfunc;
+  if(mmb1!=NULL)r.mmb1=new COperation(mmb1->Substitute(var,rop));else r.mmb1=NULL;
+  if(mmb2!=NULL)r.mmb2=new COperation(mmb2->Substitute(var,rop));else r.mmb2=NULL;
+  return r;
+}
+
+COperation COperation::Diff(const CVar& var) const //This is d / dz
+{
+  if(!ContainVar(var))return 0.0;
+  if(op==Var)return 1.0;
+  COperation **ppop1,op2;int i,j;
+  switch(op){
+  case Juxt:return(mmb1->Diff(var),mmb2->Diff(var));
+  case Add:return(mmb1->Diff(var)+mmb2->Diff(var));
+  case Sub:return(mmb1->Diff(var)-mmb2->Diff(var));
+  case Opp:return(-mmb2->Diff(var));
+  case Mult:return((*mmb1)*(mmb2->Diff(var))+(*mmb2)*(mmb1->Diff(var)));
+  case Div:if(mmb2->ContainVar(var))return(((*mmb2)*(mmb1->Diff(var))-(*mmb1)*(mmb2->Diff(var)))/((*mmb2)^2));
+  else return(mmb1->Diff(var)/(*mmb2));
+  case Pow:if(mmb2->ContainVar(var))return((*this)*(log(*mmb1)*mmb2->Diff(var)+
+						    (*mmb2)*mmb1->Diff(var)/(*mmb1)));else
+						      return (*mmb2)*mmb1->Diff(var)*((*mmb1)^(*mmb2-1));
+  case Sqrt:return(mmb2->Diff(var)/(2*sqrt(*mmb2)));
+  case NthRoot:{COperation interm=(*mmb2)^(1/(*mmb1));return interm.Diff(var);};
+  case E10:{COperation interm=(*mmb1)*(10^(*mmb2));return interm.Diff(var);};;
+  case Ln:return (mmb2->Diff(var)/(*mmb2));
+  case Exp:return (mmb2->Diff(var)*(*this));
+  case Sin:return (mmb2->Diff(var)*cos(*mmb2));
+  case Cos:return (-mmb2->Diff(var)*sin(*mmb2));
+  case Tg:return (mmb2->Diff(var)*(1+((*this)^2)));
+  case Atan:return(mmb2->Diff(var)/(1+((*mmb2)^2)));
+  case Asin:return(mmb2->Diff(var)/sqrt(1-((*mmb2)^2)));
+  case Acos:return(-mmb2->Diff(var)/sqrt(1-((*mmb2)^2)));
+  case Abs:return ((conj(*mmb2)*mmb2->Diff(var)+(*mmb2)*conj(mmb2->DiffConj(var)))/(2*abs(*mmb2)));
+  case Conj:return(conj(mmb2->DiffConj(var)));
+  case Real:return ((mmb2->Diff(var)+conj(mmb2->DiffConj(var)))/2);
+  case Imag:return ((mmb2->Diff(var)-conj(mmb2->DiffConj(var)))/double_complex(0,2));
+  case Arg:return (double_complex(0,-.5)*(mmb2->Diff(var)/(*mmb2)-conj(mmb2->DiffConj(var)/(*mmb2))));
+  case Fun:if(pfunc->type==-1||pfunc->type==0)return ErrVal;
+    if(pfunc->nvars==0)return 0.;
+    else if(pfunc->op.NMembers()>1){
+      j=pfunc->op.NMembers();
+      ppop1=new COperation*[j];
+      for(i=0;i<j;i++)ppop1[i]=new COperation(NthMember(i+1).Diff(var));
+      op2=ApplyOperator(pfunc->nvars,ppop1,&operator,);
+      for(i=0;i<pfunc->nvars;i++)delete ppop1[i];
+      delete[]ppop1;
+      return op2;  
+    }else{
+      ppop1=new COperation*[2*pfunc->nvars];
+      for(i=0;i<pfunc->nvars;i++){
+	ppop1[i]=new COperation(pfunc->op.Diff(*pfunc->ppvar[i]));
+	for(j=0;j<pfunc->nvars;j++)
+	  *ppop1[i]=ppop1[i]->Substitute(*pfunc->ppvar[j],mmb2->NthMember(j+1));
+	*ppop1[i]=(mmb2->NthMember(i+1).Diff(var))*(*ppop1[i]);
+      }
+      for(i=pfunc->nvars;i<2*pfunc->nvars;i++){
+	  ppop1[i]=new COperation(pfunc->op.DiffConj(*pfunc->ppvar[i-pfunc->nvars]));
+	  for(j=0;j<pfunc->nvars;j++)
+	    *ppop1[i]=ppop1[i]->Substitute(*pfunc->ppvar[j],mmb2->NthMember(j+1));
+	  *ppop1[i]=((conj(mmb2->NthMember(i+1-pfunc->nvars))).Diff(var))*(*ppop1[i]);
+      }
+      op2=ApplyOperator(2*pfunc->nvars,ppop1,&::operator+);
+      for(i=0;i<2*pfunc->nvars;i++)delete ppop1[i];
+      delete[]ppop1;
+      return op2;
+      // In the obtained expression, f' will have been replaced with its expression but f will remain pointing to itself ; this could cause some trouble if changing f afterwards
+    }
+  default:return ErrVal;
+  };
+}
+
+COperation COperation::DiffConj(const CVar& var) const // This one is d / d conj(z)
+{
+  if(!ContainVar(var))return 0.0;
+  if(op==Var)return 0.0;
+  COperation **ppop1,op2;int i,j;
+  switch(op){
+  case Juxt:return(mmb1->DiffConj(var),mmb2->DiffConj(var));
+  case Add:return(mmb1->DiffConj(var)+mmb2->DiffConj(var));
+  case Sub:return(mmb1->DiffConj(var)-mmb2->DiffConj(var));
+  case Opp:return(-mmb2->DiffConj(var));
+  case Mult:return((*mmb1)*(mmb2->DiffConj(var))+(*mmb2)*(mmb1->DiffConj(var)));
+  case Div:if(mmb2->ContainVar(var))return(((*mmb2)*(mmb1->DiffConj(var))-(*mmb1)*(mmb2->DiffConj(var)))/((*mmb2)^2));
+  else return(mmb1->DiffConj(var)/(*mmb2));
+  case Pow:if(mmb2->ContainVar(var))return((*this)*(log(*mmb1)*mmb2->DiffConj(var)+
+						    (*mmb2)*mmb1->DiffConj(var)/(*mmb1)));else
+						      return (*mmb2)*mmb1->DiffConj(var)*((*mmb1)^(*mmb2-1));
+  case Sqrt:return(mmb2->DiffConj(var)/(2*sqrt(*mmb2)));
+  case NthRoot:{COperation interm=(*mmb2)^(1/(*mmb1));return interm.DiffConj(var);};
+  case E10:{COperation interm=(*mmb1)*(10^(*mmb2));return interm.DiffConj(var);};;
+  case Ln:return (mmb2->DiffConj(var)/(*mmb2));
+  case Exp:return (mmb2->DiffConj(var)*(*this));
+  case Sin:return (mmb2->DiffConj(var)*cos(*mmb2));
+  case Cos:return (-mmb2->DiffConj(var)*sin(*mmb2));
+  case Tg:return (mmb2->DiffConj(var)*(1+((*this)^2)));
+  case Atan:return(mmb2->DiffConj(var)/(1+((*mmb2)^2)));
+  case Asin:return(mmb2->DiffConj(var)/sqrt(1-((*mmb2)^2)));
+  case Acos:return(-mmb2->DiffConj(var)/sqrt(1-((*mmb2)^2)));
+  case Abs:return ((conj(*mmb2)*mmb2->DiffConj(var)+(*mmb2)*conj(mmb2->Diff(var)))/(2*abs(*mmb2)));
+  case Conj:return(conj(mmb2->Diff(var)));
+  case Real:return ((mmb2->DiffConj(var)+conj(mmb2->Diff(var)))/2);
+  case Imag:return ((mmb2->DiffConj(var)-conj(mmb2->Diff(var)))/double_complex(0,2));
+  case Arg:return (double_complex(0,-.5)*(mmb2->DiffConj(var)/(*mmb2)-conj(mmb2->Diff(var)/(*mmb2))));
+  case Fun:if(pfunc->type==-1||pfunc->type==0)return ErrVal;
+    if(pfunc->nvars==0)return 0.;
+    else if(pfunc->op.NMembers()>1){
+      j=pfunc->op.NMembers();
+      ppop1=new COperation*[j];
+      for(i=0;i<j;i++)ppop1[i]=new COperation(NthMember(i+1).DiffConj(var));
+      op2=ApplyOperator(pfunc->nvars,ppop1,&operator,);
+      for(i=0;i<pfunc->nvars;i++)delete ppop1[i];
+      delete[]ppop1;
+      return op2;  
+    }else{
+      ppop1=new COperation*[2*pfunc->nvars];
+      for(i=0;i<pfunc->nvars;i++){
+	ppop1[i]=new COperation(pfunc->op.DiffConj(*pfunc->ppvar[i]));
+	for(j=0;j<pfunc->nvars;j++)
+	  *ppop1[i]=ppop1[i]->Substitute(*pfunc->ppvar[j],mmb2->NthMember(j+1));
+	*ppop1[i]=(conj(mmb2->NthMember(i+1).Diff(var)))*(*ppop1[i]);
+      }
+      for(i=pfunc->nvars;i<2*pfunc->nvars;i++){
+	ppop1[i]=new COperation(pfunc->op.Diff(*pfunc->ppvar[i-pfunc->nvars]));
+	for(j=0;j<pfunc->nvars;j++)
+	  *ppop1[i]=ppop1[i]->Substitute(*pfunc->ppvar[j],mmb2->NthMember(j+1));
+	*ppop1[i]=(mmb2->NthMember(i+1-pfunc->nvars).DiffConj(var))*(*ppop1[i]);
+      }
+      op2=ApplyOperator(2*pfunc->nvars,ppop1,&::operator+);
+      for(i=0;i<2*pfunc->nvars;i++)delete ppop1[i];
+      delete[]ppop1;
+      return op2;
+      // In the obtained expression, f' will have been replaced with its expression but f will remain pointing to itself ; this could cause some trouble if changing f afterwards
+    }
+  default:return ErrVal;
+  };
+}
+
+char* formatreal(double x)
+{
+  char*s=new char[20];
+  if(x==0)s[0]=0;
+  else if(x==(double)3.141592653589793238462643383279L)sprintf(s,"pi");
+  else sprintf(s,"%.8g",x);
+  return s;
+}
+
+char* formatimag(double x)
+{
+  char*s=new char[20];
+  if(x==0)s[0]=0;
+  else if(x==1)sprintf(s,"i");
+  else if(x==-1)sprintf(s,"-i");
+  else if(x==(double)3.141592653589793238462643383279L)sprintf(s,"i pi");
+  else if(x==(double)(-3.141592653589793238462643383279L))sprintf(s,"-i pi");
+  else sprintf(s,"%.8g i",x);
+  return s;
+}
+
+char* PrettyPrint(double_complex x)
+{
+  char*s=new char[30],*s1=formatreal(real(x)),*s2=formatimag(imag(x));
+  if(!strlen(s1)&&!strlen(s2))sprintf(s,"0");
+  else if(x==double_complex(0,1))sprintf(s,"i");
+  else if(!strlen(s1))sprintf(s,"(%s)",s2);
+  else if(!strlen(s2))sprintf(s,"%s",s1);
+  else if(imag(x)>0)sprintf(s,"(%s+%s)",s1,s2);
+  else sprintf(s,"(%s%s)",s1,s2);
+  delete s1, delete s2;
+  return s;
+}
+
+char* COperation::Expr() const
+{
+  char*s=NULL,*s1=NULL,*s2=NULL;int n=10;signed char f=0,g=0;
+  if(op==Fun)if(strlen(pfunc->name)>4)n+=strlen(pfunc->name)-4;
+  if(mmb1!=NULL){s1=mmb1->Expr();n+=strlen(s1);f=IsFunction(mmb1->op);}
+  if(mmb2!=NULL){s2=mmb2->Expr();n+=strlen(s2);g=IsFunction(mmb2->op);}
+  s=new char[n];
+  switch(op){
+  case Num:return PrettyPrint(ValC);
+  case Var:return CopyStr(pvar->name);
+  case Juxt:sprintf(s,"%s , %s",s1,s2);break;
+  case Add:
+    f=f||(mmb1->op==Juxt);
+    g=g||(mmb2->op==Juxt);
+    if(f&&g)sprintf(s,"(%s)+(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)+%s",s1,s2);else
+	if(g)sprintf(s,"%s+(%s)",s1,s2);else
+	  sprintf(s,"%s+%s",s1,s2);
+    break;
+  case Sub:
+    f=f||(mmb1->op==Juxt);
+    g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub);
+    if(f&&g)sprintf(s,"(%s)-(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)-%s",s1,s2);else
+	if(g)sprintf(s,"%s-(%s)",s1,s2);else
+	  sprintf(s,"%s-%s",s1,s2);
+    break;
+  case Opp:
+    if(mmb2->op==Add||mmb2->op==Sub||mmb2->op==Juxt)sprintf(s,"-(%s)",s2);else
+      sprintf(s,"-%s",s2);
+    break;
+  case Mult:
+    f=f||(mmb1->op==Juxt||mmb1->op==Add||mmb1->op==Sub||mmb1->op==Opp);
+    g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub||mmb2->op==Opp);
+    if(f&&g)sprintf(s,"(%s)*(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)*%s",s1,s2);else
+	if(g)sprintf(s,"%s*(%s)",s1,s2);else
+	  sprintf(s,"%s*%s",s1,s2);
+    break;
+  case Div:
+    f=f||(mmb1->op==Juxt||mmb1->op==Add||mmb1->op==Sub||mmb1->op==Opp||mmb1->op==Div);
+    g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub||mmb2->op==Opp||mmb2->op==Mult||mmb2->op==Div);
+    if(f&&g)sprintf(s,"(%s)/(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)/%s",s1,s2);else
+	if(g)sprintf(s,"%s/(%s)",s1,s2);else
+	  sprintf(s,"%s/%s",s1,s2);
+    break;
+  case Pow:
+    f=(mmb1->op!=Num&&mmb1->op!=Var);
+    g=(mmb2->op!=Num&&mmb2->op!=Var);
+    if(f&&g)sprintf(s,"(%s)^(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)^%s",s1,s2);else
+	if(g)sprintf(s,"%s^(%s)",s1,s2);else
+	  sprintf(s,"%s^%s",s1,s2);
+    break;
+  case Sqrt:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"sqrt(%s)",s2);
+    else sprintf(s,"sqrt %s",s2);
+    break;
+  case NthRoot:
+    f=(mmb1->op!=Num&&mmb1->op!=Var);
+    g=(mmb2->op!=Num&&mmb2->op!=Var);
+    if(f&&g)sprintf(s,"(%s)#(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)#%s",s1,s2);else
+	if(g)sprintf(s,"%s#(%s)",s1,s2);else
+	  sprintf(s,"%s#%s",s1,s2);
+    break;
+  case E10:
+    f=(mmb1->op!=Num&&mmb1->op!=Var);
+    g=(mmb2->op!=Num&&mmb2->op!=Var);
+    if(f&&g)sprintf(s,"(%s)E(%s)",s1,s2);else
+      if(f)sprintf(s,"(%s)E%s",s1,s2);else
+	if(g)sprintf(s,"%sE(%s)",s1,s2);else
+	  sprintf(s,"%sE%s",s1,s2);
+    break;
+  case Ln:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"log(%s)",s2);
+    else sprintf(s,"log %s",s2);
+    break;
+  case Exp:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"exp(%s)",s2);
+    else sprintf(s,"exp %s",s2);
+    break;
+  case Sin:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"sin(%s)",s2);
+    else sprintf(s,"sin %s",s2);
+    break;
+  case Cos:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"cos(%s)",s2);
+    else sprintf(s,"cos %s",s2);
+    break;
+  case Tg:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"tan(%s)",s2);
+    else sprintf(s,"tan %s",s2);
+    break;
+  case Atan:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"atan(%s)",s2);
+    else sprintf(s,"atan %s",s2);
+    break;
+  case Asin:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"asin(%s)",s2);
+    else sprintf(s,"asin %s",s2);
+    break;
+  case Acos:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"acos(%s)",s2);
+    else sprintf(s,"acos %s",s2);
+    break;
+  case Abs:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"abs(%s)",s2);
+    else sprintf(s,"abs %s",s2);
+    break;
+  case Real:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"real(%s)",s2);
+    else sprintf(s,"real %s",s2);
+    break;
+  case Imag:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"imag(%s)",s2);
+    else sprintf(s,"imag %s",s2);
+    break;
+  case Conj:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"conj(%s)",s2);
+    else sprintf(s,"conj %s",s2);
+    break;
+  case Arg:
+    g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+    if(g)sprintf(s,"arg(%s)",s2);
+    else sprintf(s,"arg %s",s2);
+    break;
+  case Fun:
+    sprintf(s,"%s(%s)",pfunc->name,s2);
+    break;
+  default:return CopyStr("Error");
+  };
+  if(s1!=NULL)delete[] s1;if(s2!=NULL)delete[] s2;
+  return s;
+}
+
+const double sqrtmaxfloat=sqrt(DBL_MAX);
+const double sqrtminfloat=sqrt(DBL_MIN);
+
+void Addition(double_complex*&p)
+{if(*p==ErrVal||norm(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||norm(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p+=(*(p+1));}
+void Soustraction(double_complex*&p)
+{if(*p==ErrVal||norm(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||norm(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p-=(*(p+1));}
+void Multiplication(double_complex*&p)
+{if(norm(*p)<sqrtminfloat){*--p=0;return;};
+if(*p==ErrVal||norm(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(norm(*(--p))<sqrtminfloat){*p=0;return;};
+if(*p==ErrVal||norm(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p*=(*(p+1));}
+void Division(double_complex*&p)
+{if(norm(*p)<sqrtminfloat||*p==ErrVal||norm(*p)>sqrtmaxfloat)
+  {*(--p)=ErrVal;return;};
+if(norm(*(--p))<sqrtminfloat)*p=0;else if(*p==ErrVal||norm(*p)>sqrtmaxfloat)
+  {*p=ErrVal;return;};
+ *p/=(*(p+1));}
+void Puissance(double_complex*&p)
+{double_complex v2=*p--,v1=*p;
+if(v2==ErrVal||v1==ErrVal||norm(v2)*log(norm(v1))>DBL_MAX_EXP){*p=ErrVal;return;};
+*p=(v1==0.?0:pow(v1,v2));}
+void RacineN(double_complex*&p)
+{double_complex v2=*p--,v1=*p;
+if(v1==ErrVal||v2==ErrVal||v1==0.||norm(v2)*log(norm(v1))<DBL_MIN_EXP){*p=ErrVal;return;};
+*p=pow(v2,double_complex(1,0)/v1);return;}
+void Puiss10(double_complex*&p)
+{if(norm(*p)<sqrtminfloat){*(--p)=0;return;};
+if(*p==ErrVal||norm(*p)>DBL_MAX_10_EXP){*(--p)=ErrVal;return;};
+if(norm(*(--p))<sqrtminfloat)*p=0;else if(*p==ErrVal||norm(*p)>sqrtmaxfloat)
+	{*p=ErrVal;return;};
+*p*=pow(double_complex(10.),*(p+1));}
+void NextVal(double_complex*&){}
+void  RFunc(double_complex*&){}
+void  JuxtF(double_complex*&){}
+void Absolu(double_complex*&p){*p=((*p==ErrVal)?ErrVal:norm(*p));}
+void Oppose(double_complex*&p){*p=((*p==ErrVal)?ErrVal:-*p);}
+void Reelle(double_complex*&p){*p=((*p==ErrVal)?ErrVal:real(*p));}
+void Imaginaire(double_complex*&p){*p=((*p==ErrVal)?ErrVal:imag(*p));}
+void Conjugue(double_complex*&p){*p=((*p==ErrVal)?ErrVal:conj(*p));}
+void ArcSinus(double_complex*&p)
+{*p=((*p==ErrVal)?ErrVal:asin(*p));}
+void ArcCosinus(double_complex*&p)
+{*p=((*p==ErrVal||*p==double_complex(0,1)||*p==double_complex(0,-1))?ErrVal:acos(*p));}
+void ArcTangente(double_complex*&p)
+{*p=((*p==ErrVal)?ErrVal:atan(*p));}
+void Logarithme(double_complex*&p)
+{*p=((*p==ErrVal||*p==0.)?ErrVal:log(*p));}
+void Argument(double_complex*&p)
+{*p=((*p==ErrVal||*p==0.)?ErrVal:arg(*p));}
+void Exponentielle(double_complex*&p)
+{*p=((*p==ErrVal||norm(*p)>DBL_MAX_EXP)?ErrVal:exp(*p));}
+void Sinus(double_complex*&p)
+{*p=((*p==ErrVal||norm(*p)>DBL_MAX_EXP)?ErrVal:sin(*p));}
+void Tangente(double_complex*&p)
+{*p=((*p==ErrVal||norm(*p)>DBL_MAX_EXP)?ErrVal:tan(*p));}
+void Cosinus(double_complex*&p)
+{*p=((*p==ErrVal||norm(*p)>DBL_MAX_EXP)?ErrVal:cos(*p));}
+void Racine(double_complex*&p)
+{*p=((*p==ErrVal||norm(*p)>sqrtmaxfloat)?ErrVal:sqrt(*p));}
+void FonctionError(double_complex*&p){*p=ErrVal;}
+inline void ApplyRFunc(PCFunction rf,double_complex*&p)
+{p-=rf->nvars-1;*p=rf->Val(p);}
+
+double_complex COperation::Val() const
+{
+  pfoncld*p1=pinstr;double_complex**p2=pvals,*p3=ppile-1;PCFunction*p4=pfuncpile;
+  for(;*p1!=NULL;p1++)
+    if(*p1==&NextVal)*(++p3)=**(p2++);else 
+      if(*p1==&RFunc) ApplyRFunc(*(p4++),p3);
+      else (**p1)(p3);
+  return *p3;
+}
+
+void BCDouble(pfoncld*&pf,pfoncld*pf1,pfoncld*pf2,
+	      double_complex**&pv,double_complex**pv1,double_complex**pv2,
+	      double_complex*&pp,double_complex*pp1,double_complex*pp2,
+	      CFunction**&prf,CFunction**prf1,CFunction**prf2,
+	      pfoncld f)
+{
+  pfoncld*pf3,*pf4=pf1;long n1,n2;
+  for(n1=0;*pf4!=NULL;pf4++,n1++);for(n2=0,pf4=pf2;*pf4!=NULL;pf4++,n2++);
+  pf=new pfoncld[n1+n2+2];
+  for(pf3=pf,pf4=pf1;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  for(pf4=pf2;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  *pf3++=f;*pf3=NULL;//delete[]pf1,pf2;
+  double_complex**pv3,**pv4=pv1;
+  for(n1=0;*pv4!=NULL;pv4++,n1++);for(n2=0,pv4=pv2;*pv4!=NULL;pv4++,n2++);
+  pv=new double_complex*[n1+n2+1];
+  for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  for(pv4=pv2;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  *pv3=NULL;//delete[]pv1,pv2;
+  double_complex*pp3,*pp4=pp1;
+  for(n1=0;*pp4!=ErrVal;pp4++,n1++);for(n2=0,pp4=pp2;*pp4!=ErrVal;pp4++,n2++);
+  pp=new double_complex[n1+n2+1];  // Really need to add and not to take max(n1,n2) in case of Juxt operator
+  for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  for(pp4=pp2;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  *pp3=ErrVal;//delete[]pp1,pp2;
+  PCFunction*prf3,*prf4=prf1;
+  for(n1=0;*prf4!=NULL;prf4++,n1++);for(n2=0,prf4=prf2;*prf4!=NULL;prf4++,n2++);
+  prf=new PCFunction[n1+n2+1];
+  for(prf3=prf,prf4=prf1;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  for(prf4=prf2;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  *prf3=NULL;//delete[]prf1,prf2;
+}
+
+void BCSimple(pfoncld*&pf,pfoncld*pf1,double_complex**&pv,double_complex**pv1,
+	double_complex*&pp,double_complex*pp1,CFunction**&prf,CFunction**prf1,pfoncld f)
+{
+  pfoncld*pf3,*pf4=pf1;long n;
+  for(n=0;*pf4!=NULL;pf4++,n++);
+  pf=new pfoncld[n+2];
+  for(pf4=pf1,pf3=pf;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  *pf3++=f;*pf3=NULL;//delete[]pf1;
+  double_complex**pv3,**pv4=pv1;
+  for(n=0;*pv4!=NULL;pv4++,n++);
+  pv=new double_complex*[n+1];
+  for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  *pv3=NULL;//delete[]pv1;
+  double_complex*pp3,*pp4=pp1;
+  for(n=0;*pp4!=ErrVal;pp4++,n++);
+  pp=new double_complex[n+1];
+  for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  *pp3=ErrVal;//delete[]pp1;
+  CFunction**prf3,**prf4=prf1;
+  for(n=0;*prf4!=NULL;prf4++,n++);
+  prf=new CFunction*[n+1];
+  for(prf3=prf,prf4=prf1;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  *prf3=NULL;//delete[]prf1;
+}
+
+void BCFun(pfoncld*&pf,pfoncld*pf1,double_complex**&pv,double_complex**pv1,
+	double_complex*&pp,double_complex*pp1,CFunction**&prf,CFunction**prf1,PCFunction rf)
+{
+  pfoncld*pf3,*pf4=pf1;long n;
+  for(n=0;*pf4!=NULL;pf4++,n++);
+  pf=new pfoncld[n+2];
+  for(pf4=pf1,pf3=pf;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+  *pf3++=&RFunc;*pf3=NULL;//delete[]pf1;
+  double_complex**pv3,**pv4=pv1;
+  for(n=0;*pv4!=NULL;pv4++,n++);
+  pv=new double_complex*[n+1];
+  for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+  *pv3=NULL;//delete[]pv1;
+  double_complex*pp3,*pp4=pp1;
+  for(n=0;*pp4!=ErrVal;pp4++,n++);
+  pp=new double_complex[n+1];
+  for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+  *pp3=ErrVal;//delete[]pp1;
+  PCFunction*prf3,*prf4=prf1;
+  for(n=0;*prf4!=NULL;prf4++,n++);
+  prf=new PCFunction[n+2];
+  for(prf4=prf1,prf3=prf;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+  *prf3++=rf;*prf3=NULL;//delete[]pf1;
+}
+
+void COperation::BuildCode()
+{
+  //  if(mmb1!=NULL)mmb1->BuildCode();if(mmb2!=NULL)mmb2->BuildCode();
+  if(pinstr!=NULL){delete[]pinstr;pinstr=NULL;}
+  if(pvals!=NULL){delete[]pvals;pvals=NULL;}
+  if(ppile!=NULL){delete[]ppile;ppile=NULL;}
+  if(pfuncpile!=NULL){delete[]pfuncpile;pfuncpile=NULL;}
+  switch(op){
+  case ErrOp:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+    pvals=new double_complex*[2];pvals[0]=new double_complex(ErrVal);pvals[1]=NULL;
+    ppile=new double_complex[2];ppile[0]=0;ppile[1]=ErrVal;
+    pfuncpile=new CFunction*[1];pfuncpile[0]=NULL;
+    break;
+  case Num:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+    pvals=new double_complex*[2];pvals[0]=new double_complex(ValC);pvals[1]=NULL;
+    ppile=new double_complex[2];ppile[0]=0;ppile[1]=ErrVal;
+    pfuncpile=new CFunction*[1];pfuncpile[0]=NULL;break;
+  case Var:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+    pvals=new double_complex*[2];pvals[0]=pvarval;pvals[1]=NULL;
+    ppile=new double_complex[2];ppile[0]=0;ppile[1]=ErrVal;
+    pfuncpile=new CFunction*[1];pfuncpile[0]=NULL;break;
+  case Juxt:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+		     pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&JuxtF);
+  break;case Add:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+		     pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Addition);
+  break;case Sub:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Soustraction);
+  break;case Mult:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			   pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Multiplication);
+  break;case Div:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Division);
+  break;case Pow:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Puissance);
+  break;case NthRoot:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			      pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&RacineN);
+  break;case E10:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+			  pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Puiss10);
+  break;case Opp:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Oppose);
+  break;case Sin:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Sinus);
+  break;case Sqrt:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Racine);
+  break;case Ln:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			 ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Logarithme);
+  break;case Exp:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Exponentielle);
+  break;case Cos:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Cosinus);
+  break;case Tg:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			 ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Tangente);
+  break;case Atan:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&ArcTangente);
+  break;case Asin:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&ArcSinus);
+  break;case Acos:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&ArcCosinus);
+  break;case Abs:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Absolu);
+  break;case Real:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Reelle);
+  break;case Imag:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Imaginaire);
+  break;case Conj:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			   ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Conjugue);
+  break;case Arg:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			  ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Argument);
+  break;case Fun:BCFun(pinstr,mmb2->pinstr,pvals,mmb2->pvals,ppile,
+		       mmb2->ppile,pfuncpile,mmb2->pfuncpile,pfunc);
+  break;default:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+			 ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&FonctionError);
+  }
+}
diff --git a/mathexpr/mathexpr/mathexpr_c.h b/mathexpr/mathexpr/mathexpr_c.h
new file mode 100644
index 0000000..faa4970
--- /dev/null
+++ b/mathexpr/mathexpr/mathexpr_c.h
@@ -0,0 +1,160 @@
+/*
+
+mathexpr_c.h version 2.0
+
+Copyright (c) 1997-2000 Yann OLLIVIER
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+*/
+
+#ifndef __MATHEXPR_C_H
+#define __MATHEXPR_C_H
+
+#include<string.h>
+#include<stdio.h>
+#include<stdlib.h>
+#include<complex>
+#include<float.h>
+
+#ifndef __STD_COMPLEX_
+#define double_complex std::complex<double>
+double_complex tan(double_complex);
+double_complex acos(double_complex);
+double_complex asin(double_complex);
+double_complex atan(double_complex);
+#endif
+
+// Compatibility with long double-typed functions
+#define atanl atan
+#define asinl asin
+#define acosl acos
+#define expl exp
+#define logl log
+#define powl pow
+#define pow10l(x) pow(10,x)
+#define fabsl fabs
+#define cosl cos
+#define sinl sin
+#define tanl tan
+#define fmodl fmod
+#define sqrtl sqrt
+
+
+// Warning : everything will fail with strings longer than 32767 chars
+
+const double_complex ErrVal=double_complex(DBL_MAX,0);
+
+//Class definitions for operations
+
+class CVar{
+ public:
+  char*name;double_complex*pval;
+  CVar(const CVar&);
+  CVar(const char*,double_complex*);
+  ~CVar();
+  friend int operator==(const CVar&,const CVar&);
+};
+
+typedef CVar* PCVar;
+
+enum COperator{ErrOp,Juxt,Num,Var,Add,Sub,Opp,Mult,Div,Pow,Sqrt,
+	       NthRoot,Real,Imag,Conj,Abs,Arg,Sin,Cos,Tg,Ln,Exp,Acos,Asin,Atan,E10,Fun};
+
+typedef void ((*pfoncld)(double_complex*&));
+
+class COperation;
+typedef COperation* PCOperation;
+class CFunction;
+typedef CFunction* PCFunction;
+
+class COperation{
+  pfoncld*pinstr;double_complex**pvals;double_complex*ppile;CFunction**pfuncpile;
+  mutable signed char containfuncflag;
+  void BuildCode();
+	void Destroy();
+ public:
+  COperator op;
+  PCOperation mmb1,mmb2;
+  double_complex ValC;const CVar* pvar;double_complex*pvarval;
+  CFunction* pfunc;
+  COperation();
+  COperation(const COperation&);
+  COperation(double);COperation(double_complex);
+  COperation(const CVar&);
+  COperation(char*sp,int nvarp=0,PCVar*ppvarp=NULL,int nfuncp=0,PCFunction*ppfuncp=NULL);
+  ~COperation();
+  double_complex Val() const;
+  signed char ContainVar(const CVar&) const;
+  signed char ContainFunc(const CFunction&) const;
+  signed char ContainFuncNoRec(const CFunction&) const; // No recursive test on subfunctions
+  COperation NthMember(int) const;int NMembers() const;
+  signed char HasError(const COperation* =NULL) const;
+  COperation& operator=(const COperation&);
+  friend int operator==(const COperation&,const double);
+  friend int operator==(const COperation&,const double_complex);
+  friend int operator==(const COperation&,const COperation&);
+  friend int operator!=(const COperation&,const COperation&);
+  COperation operator+() const;COperation operator-() const;
+  friend COperation operator,(const COperation&,const COperation&);
+  friend COperation operator+(const COperation&,const COperation&);
+  friend COperation operator-(const COperation&,const COperation&);
+  friend COperation operator*(const COperation&,const COperation&);
+  friend COperation operator/(const COperation&,const COperation&);
+  friend COperation operator^(const COperation&,const COperation&);  // Caution: wrong associativity and precedence
+  friend COperation real(const COperation&);
+  friend COperation imag(const COperation&);
+  friend COperation conj(const COperation&);
+  friend COperation arg(const COperation&);
+  friend COperation sqrt(const COperation&);
+  friend COperation abs(const COperation&);
+  friend COperation sin(const COperation&);
+  friend COperation cos(const COperation&);
+  friend COperation tan(const COperation&);
+  friend COperation log(const COperation&);
+  friend COperation exp(const COperation&);
+  friend COperation acos(const COperation&);
+  friend COperation asin(const COperation&);
+  friend COperation atan(const COperation&);
+  friend COperation ApplyOperator(int,COperation**,COperation (*)(const COperation&,const COperation&));
+  COperation Diff(const CVar&) const; //  Differentiate w.r.t a variable
+  COperation DiffConj(const CVar&) const; // This one is d / d conj(z)
+  char* Expr() const;
+  COperation Substitute(const CVar&,const COperation&) const;
+};
+
+class CFunction{
+  double_complex*buf;
+public:
+  signed char type;
+  double_complex ((*pfuncval)(double_complex));
+  COperation op;int nvars;CVar** ppvar;
+  char*name;
+  CFunction();
+  CFunction(double_complex ((*pfuncvalp)(double_complex)));
+  CFunction(const COperation& opp, CVar* pvarp);
+  CFunction(const COperation& opp, int nvarsp, CVar**ppvarp);
+  CFunction(const CFunction&);
+  ~CFunction();
+  CFunction& operator=(const CFunction&);
+  void SetName(const char*s);
+  double_complex Val(double_complex) const;
+  double_complex Val(double_complex*) const;
+  friend int operator==(const CFunction&,const CFunction&);
+  COperation operator()(const COperation&);
+};
+
+char* PrettyPrint(double_complex);
+
+char* MidStr(const char*s,int i1,int i2);
+char* CopyStr(const char*s);
+char* InsStr(const char*s,int n,const char c);
+signed char EqStr(const char*s,const char*s2);
+signed char CompStr(const char*s,int n,const char*s2);
+char* DelStr(const char*s,int n);
+
+#endif