{
    $id:                                                       $
    Copyright (c) 2000 by Marco van de Voort (marco@freepascal.org)
    member of the Free Pascal development team

    See the file COPYING.FPC, included in this distribution,
    for details about the copyright. (LGPL)

    TExpression class which does symbolic manipulation.

    Derivation routine based on 20 lines of conceptual pseudo code
    provided by Osmo Ronkanen.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 **********************************************************************

Problems:

- Often untested. I used my HP48g to check some complex derivatives,
   but more thorough checking and errorhandling is necessary.
- RPN to Infix adds ()'s when not necessary. Should be made aware of precedence.
 (partially fixed)
}

{Should be moved to a math unit. Calculate x! with a 23 degree polynomal for
 ArbFloat values. From Numlib.


Works for extended only. Redefine TC1 and TC3 for your numeric type
if you want to use something else.}

type Float10Arb =ARRAY[0..9] OF BYTE;

const

   TC1 : Float10Arb = (0,0,$00,$00,$00,$00,0,128,192,63);         {Eps}
   TC3 : Float10Arb = (1,0,0,0,0,0,0,0,0,0);                      {3.64519953188247460E-4951}


var {Looks ugly, but is quite handy.}
    macheps  : ArbFloat absolute TC1;  { macheps = r - 1,  with r
                                        the smallest ArbFloat > 1}
    midget   : ArbFloat absolute TC3;  { the smallest positive ArbFloat}


function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
var   pa : ^ArbFloat; {FPC extension. Uses ^ some array of ArbFloat in TP}
       i : ArbInt;
    polx : ArbFloat;
begin
  pa:=@a;
  polx:=0;
  for i:=n downto 0 do
    polx:=polx*x+pa[i]; {and pa^[i] here}
  spepol:=polx
end {spepol};


function spegam(x: ArbFloat): ArbFloat;
const

    tmax = 170;
    a: array[0..23] of ArbFloat =
    ( 8.86226925452758013e-1,  1.61691987244425092e-2,
      1.03703363422075456e-1, -1.34118505705967765e-2,
      9.04033494028101968e-3, -2.42259538436268176e-3,
      9.15785997288933120e-4, -2.96890121633200000e-4,
      1.00928148823365120e-4, -3.36375833240268800e-5,
      1.12524642975590400e-5, -3.75499034136576000e-6,
      1.25281466396672000e-6, -4.17808776355840000e-7,
      1.39383522590720000e-7, -4.64774927155200000e-8,
      1.53835215257600000e-8, -5.11961333760000000e-9,
      1.82243164160000000e-9, -6.13513953280000000e-10,
      1.27679856640000000e-10,-4.01499750400000000e-11,
      4.26560716800000000e-11,-1.46381209600000000e-11);

var tvsmall, t, g: ArbFloat;
             m, i: ArbInt;
begin
  if sizeof(ArbFloat) = 6
  then
    tvsmall:=2*midget
  else
    tvsmall:=midget;
  t:=abs(x);
  if t > tmax
  then
    RunError(407);
  if t < macheps
  then
    begin
      if t < tvsmall
      then
        RunError(407);
      spegam:=1/x
    end
  else  { abs(x) >= macheps }
    begin
      m:=trunc(x);
      if x > 0
      then
        begin
          t:=x-m; m:=m-1; g:=1;
          if m<0
          then
            g:=g/x
          else
            if m>0
            then
              for i:=1 to m do
                g:=(x-i)*g
        end
      else { x < 0 }
        begin
          t:=x-m+1;
          if t=1
          then
            RunError(407);
          m:=1-m;
          g:=x;
          for i:=1 to m do
            g:=(i+x)*g;
          g:=1/g
        end;
      spegam:=spepol(2*t-1, a[0], sizeof(a) div sizeof(ArbFloat) - 1)*g
    end { abs(x) >= macheps }
end; {spegam}


Procedure ExprInternalError(A,B:ArbInt);

VAR S,S2 : String;

begin
 Str(ORD(A),S);
 Str(ORD(B),S2);
 Raise EExprIE.Create(SExprIE+S+' '+S2);
end;

CONSTRUCTOR TExpression.Create(Infix:String);

var dummy:String;

begin
 ExprTree:=NIL;

 if (Infix<>'') then
  ExprTree:=InfixToParseTree(Infix,Dummy);
 InfixCache:=Infix;
 InfixClean:=True;                      {Current pnode status' infix
                                            representation is in infixcache}
end;

CONSTRUCTOR TExpression.EmptyCreate;

begin
 ExprTree:=Nil;
 InfixClean:=false;
end;

Procedure TExpression.SetNewInfix(Infix:String);

var dummy:String;

begin
 if Assigned(ExprTree) Then
  Dispose(ExprTree);
 if infix<>'' then
  ExprTree:=InFixToParseTree(Infix,Dummy)
 else
  ExprTree:=NIL;
 InfixClean:=True;
 InfixCache:=Infix;
end;

Destructor TExpression.Destroy;

begin
 If assigned(ExprTree) then
  DisposeExpr(ExprTree);
 inherited Destroy;
end;

function TExpression.GetRPN :String;

begin
 if ExprTree=NIL Then
  Result:='0'
 else
  Result:=ParseTreeToRpn(ExprTree);
end;

function TExpression.GetInfix:String;

begin
 if Not InfixClean then
  begin
   If ExprTree=NIL THEN
    InfixCache:='0'
   else
    InfixCache:=ParseTreeToInfix(ExprTree);
   InfixClean:=True;
  end;
 Result:=InfixCache;
end;

Function  TExpression.GetIntValue:LongInt;

begin
 SimplifyConstants;
 If ExprTree^.NodeType<>Iconstnode then
  Raise ENotInt.Create(SExprNotInt);
 result:=ExprTree^.ivalue;
end;

Procedure TExpression.SetIntValue(val:Longint);

begin
 if ExprTree<> NIL then
  DisposeExpr(ExprTree);
 New(ExprTree);
 ExprTree^.NodeType:=iconstnode;
 ExprTree^.Ivalue:=Val;
end;

Function  TExpression.GetFloatValue:ArbFloat;

begin
 If ExprTree^.NodeType<>constnode then
  Raise ENotFloat.Create(SExprNotFloat);
 result:=ExprTree^.value;
end;

Procedure TExpression.SetFloatValue(val:ArbFloat);

begin
 if ExprTree<> NIL then
  DisposeExpr(ExprTree);
 New(ExprTree);
 ExprTree^.NodeType:=constnode;
 ExprTree^.value:=Val;
end;

procedure TExpression.Simpleop(expr:TExpression;oper:calcop);

begin
 exprtree:=NewCalc(oper,exprtree,CopyTree(expr.exprtree));
 InFixCache:='garbadge';
 InfixClean:=False;
end;

function TExpression.Simpleopwithresult(expr:TExpression;oper:calcop):TExpression;

var tmp:pnode;

begin
 result.EmptyCreate;
 result.SimplificationLevel:=simplificationlevel;
 result.exprtree:=NewCalc(Oper,CopyTree(ExprTree),CopyTree(Expr.ExprTree));
end;

procedure TExpression.Addto(Expr:TExpression);

begin
 simpleop(expr,addo);
end;


procedure TExpression.SubFrom(Expr:TExpression);

begin
 simpleop(expr,subo);
end;

procedure TExpression.Times(Expr:texpression);

begin
 simpleop(expr,mulo);
end;

procedure TExpression.Divby(Expr:TExpression);

begin
 simpleop(expr,dvdo);
end;


procedure TExpression.RaiseTo(Expr:TExpression);

begin
 simpleop(expr,powo);
end;


function TExpression.add(Expr:TExpression):TExpression;


begin
 result:=Simpleopwithresult(expr,addo);
end;

function TExpression.sub(Expr:TExpression):TExpression;


begin
 result:=Simpleopwithresult(expr,subo);
end;

function TExpression.dvd(Expr:TExpression):TExpression;

begin
 result:=Simpleopwithresult(expr,dvdo);
end;


function TExpression.mul(Expr:TExpression):TExpression;

begin
 result:=Simpleopwithresult(expr,mulo);
end;

Function  TExpression.IntDerive(const derivvariable:String;theexpr:pnode):pnode;

function Deriv(t:pnode):pnode;
{Derive subexpression T. Returns NIL if subexpression derives to 0, to avoid
unnecessary (de)allocations. This is the reason why NewCalc is so big.}

var x     : ArbFloat;
    p1,p2 : pnode;

begin
  Deriv:=nil;
  if (t=nil) then               {Early out}
   exit;
  with t^ do begin
   case nodetype of
     VarNode:    if upcase(variable)=derivvariable then
                  Deriv:=NewiConst(ArbInt(1))
                 else
                  Deriv:=NIL;
     ConstNode : Deriv:=NIL;
     IConstNode: Deriv:=NIL;
     CalcNode:  begin
                  case op of
                   addo,
                   subo:  Deriv:=NewCalc(op,Deriv(left),Deriv(right));
                   mulo:  Deriv:=NewCalc(addo,
                                 NewCalc(mulo,Deriv(left),copyTree(right)),
                                 NewCalc(mulo,Deriv(right),copytree(left)));

                   dvdo:  Deriv:=NewCalc(dvdo,
                                  NewCalc(subo,
                                   NewCalc(mulo,Deriv(left),copyTree(right)),
                                   NewCalc(mulo,Deriv(right),copytree(left))),
                                   NewFunc(sqrx,CopyTree(right)));
                  powo: begin
                         p1:=Deriv(Right);
                         if P1<>NIL then
                          p1:=NewCalc(mulo,p1,NewFunc(Lnx,CopyTree(Left))); { ln(l)*r'}
                         p2:=Deriv(Left);
                         if P2<>NIL then
                          p2:=Newcalc(Mulo,CopyTree(Right),newcalc(mulo,p2,
                                newfunc(invx,CopyTree(left))));
                         IF (P1<>nil) and (p2<>nil) then
                          deriv:=newcalc(mulo,CopyTree(t),newcalc(addo,p1,p2))
                         else
                          if (P1=NIL) and (P2=NIL) then {Simplify first to avoid this!}
                           deriv:=NIL
                          else
                           begin
                            if P1=NIL THEN { one of both is constant}
                             P1:=P2;       {The appopriate term is now in P1}
                            deriv:=newcalc(mulo,CopyTree(t),p1);
                           end;
                        end;
                     end;
                  end;
     FuncNode: begin
                 case fun of
                   invx:  Deriv:=NewCalc(dvdo,
					  NewFunc(Minus,Deriv(son)),
		          		  NewFunc(sqrx,CopyTree(son)));

                   minus: Deriv:=NewFunc(minus,Deriv(son));
                   sinx:  Deriv:=NewCalc(Mulo,
                            NewFunc(cosx,Copytree(son)),
                            Deriv(son));
                   cosx:  deriv:=NewCalc(mulo,
                            NewFunc(minus,NewFunc(sinx,copytree(son))),
                            Deriv(son));
                   tanx:  deriv:=Newcalc(dvdo,deriv(son),
                            newfunc(sqrx,newfunc(cosx,copytree(son))));
                   sqrx:  deriv:=newcalc(mulo, newiconst(2),
                                  newcalc(mulo,copytree(son),deriv(son)));
                                  { dx*1 /(2*sqrt(x)) }
                   sqrtx: deriv:=newcalc(mulo, deriv(son),newcalc(dvdo,newiconst(1),
                                       newcalc(mulo,newiconst(2),newfunc(sqrtx,copytree(son)))));
                   lnx :  deriv:=newcalc(mulo,newcalc(dvdo,newiconst(1),CopyTree(son)),
                                             deriv(son)); { dln(x)=x' * 1/x}
                   expx:  deriv:=newcalc(mulo,newfunc(expx,copytree(son)),deriv(son));
                 cotanx:  deriv:=newfunc(minus,Newcalc(dvdo,deriv(son),  { -dx/sqr(sin(x))}
                            newfunc(sqrx,newfunc(sinx,copytree(son)))));
                  coshx:  deriv:=newcalc(mulo,newfunc(sinhx,copytree(son)),deriv(son));
                  sinhx:  deriv:=newcalc(mulo,newfunc(coshx,copytree(son)),deriv(son));
               arcsinhx, {according to HP48?}
                arcsinx:  deriv:=newcalc(dvdo,deriv(son),newfunc(sqrtx,newcalc(subo,newiconst(1),
                                newfunc(sqrx,copytree(son)))));
                arccosx:  deriv:=newfunc(minus,newcalc(dvdo,deriv(son),
                                        newfunc(sqrtx,newcalc(subo,newiconst(1),newfunc(sqrx,copytree(son))))));
                arctanx:  deriv:=newcalc(dvdo,deriv(son),newcalc(addo,newiconst(1),newfunc(sqrx,copytree(son))));
                 log10x:  deriv:=newcalc(mulo,newcalc(dvdo,newconst(0.434294481902),CopyTree(son)),
                                             deriv(son)); { dlog10(x)=x' * log10(e)/x}
                  log2x:  deriv:=newcalc(mulo,newcalc(dvdo,newconst(1.44269504089),CopyTree(son)),
                                             deriv(son)); { dlog2(x)=x' * log2(e)/x}
                  stepx:   ;  {Should raise exception, not easily derivatable}
                  tanhx:  deriv:=newcalc(dvdo,deriv(son),newfunc(sqrx,newfunc(coshx,copytree(son))));
               arctanhx:  deriv:=newcalc(dvdo,deriv(son),newfunc(sqrtx,newcalc(addo,newiconst(1),
                                           newfunc(sqrx,copytree(son)))));
               arccoshx:  deriv:=NewCalc(dvdo,deriv(son),newcalc(mulo,newcalc(subo,newfunc(sqrtx,copytree(son)),newiconst(1)),
                                                                    newcalc(addo,newfunc(sqrtx,copytree(son)),newiconst(1))));
     lnxpix,arctan2x,
     hypotx,lognx : ; {Should also raise exceptions, not implemented yet}
                    end;
                  end;
     Func2Node: begin
                  if son2left^.nodetype=constnode then
                   x:=son2left^.value
                  else
                   x:=son2left^.ivalue;
                  case fun of
                 lognx:  deriv:=newcalc(mulo,newcalc(dvdo,newconst(logn(x,2.71828182846)),
                                 CopyTree(son2right)),deriv(son2right));
                                 { dlogn(x)=x' * log(n,e)/x}
                 Powerx: begin
                         p1:=Deriv(Son2Right);
                         if P1<>NIL then
                          p1:=NewCalc(mulo,p1,NewFunc(Lnx,CopyTree(Son2Left))); { ln(l)*r'}
                         p2:=Deriv(Son2Left);
                         if P2<>NIL then
                          p2:=Newcalc(Mulo,CopyTree(Son2Right),newcalc(mulo,p2,
                                newfunc(invx,CopyTree(Son2Left))));
                         IF (P1<>nil) and (p2<>nil) then
                          deriv:=newcalc(mulo,CopyTree(t),newcalc(addo,p1,p2))
                         else
                          if (P1=NIL) and (P2=NIL) then {Simplify first to avoid this!}
                           deriv:=NIL
                          else
                           begin
                            if P1=NIL THEN { one of both is constant}
                             P1:=P2;       {The appopriate term is now in P1}
                            deriv:=newcalc(mulo,CopyTree(t),p1);
                           end;
                        end;
                     end;
                end;
           end;
    end; {WITH}
end;

begin
 Result:=Deriv(theexpr);
end;

function TExpression.power(Expr:TExpression):TExpression;


begin
 result:=Simpleopwithresult(expr,powo);
end;


Function  TExpression.Derive(derivvariable:String):TExpression;

var tmpvar : Pnode;
    DerivObj: TExpression;

begin
 derivvariable:=upcase(derivvariable);
 Tmpvar:=intDerive(derivvariable,exprtree);

 DerivObj:=TExpression.emptycreate;
 If tmpvar=NIL then
  derivobj.ExprTree:=NewIconst(0)
 else
  derivobj.exprtree:=tmpvar;
 derivobj.simplificationlevel:=simplificationlevel;
 DerivObj.InfixClean:=False;
 result:=derivobj;
end;

function ipower(x,y:ArbInt):ArbInt;

var tmpval : ArbInt;

begin
 if y<0 then
  ; {exception}
 if y=0 then
   result:=1
 else
  begin
   result:=x;
   if y<>1 then
    for tmpval:=2 to y do
     result:=result*x;
  end;
end;

function ifacul(x:ArbInt):ArbInt;

var tmpval : ArbInt;

begin
 if x<0 then
  ; {exception}
 if x=0 then
   result:=1
 else
  begin
   result:=1;
   if x<>1 then
    for tmpval:=2 to x do
     result:=result*tmpval;
  end;
end;

function EvaluateFunction(funcname:funcop;param:ArbFloat):ArbFloat;

var Intermed : integer;

begin
      case funcname of
       cosx : result:=Cos(param);
       sinx : result:=sin(param);
       tanx : result:=tan(param);
       sqrx : result:=sqr(param);
      sqrtx : result:=sqrt(param);
       expx : result:=exp(param);
        lnx : result:=ln(param);
     cotanx : result:=cotan(param);
    arcsinx : result:=arcsin(param);
    arccosx : result:=arccos(param);
    arctanx : result:=arctan(param);
      sinhx : result:=sinh(param);
      coshx : result:=cosh(param);
      tanhx : result:=tanh(param);
   arcsinhx : result:=arcsinh(param);
   arccoshx : result:=arccosh(param);
   arctanhx : result:=arctanh(param);
     log10x : result:=log10(param);
      log2x : result:=log2(param);
     lnxpix : result:=lnxp1(param);
     faculx : result:=spegam(param+1.0);
         else
          ExprInternalError(2,ord(funcname));
      end;
   If Result<1E-4900 then {Uncertainty in sinus(0.0)}
    Result:=0;
end;


procedure TExpression.SimplifyConstants;

//procedure internalsimplify (expr:pnode;InCalc:Boolean;parent:pnode);
procedure internalsimplify (expr:pnode);

function isconst(p:pnode):boolean;

begin
 isconst:=(p^.nodetype=iconstnode) or (p^.nodetype=constnode);
end;

function isconstnil(p:pnode):boolean;
begin
 IsConstNil:=false;
 if (p^.nodetype=iconstnode) and (P^.ivalue=0) then
  IsConstNil:=True;
 If (p^.nodetype=constnode) and (P^.value=0) then
  IsConstNil:=True
end;

var val1,val2 : ArbFloat;
    ival1,
    ival2 : ArbInt;

function setupoperation(operat:calcop;simlevel:longint;Postprocess:boolean;param2func:boolean):longint;

function dosimple(mode:longint;theleft,theright:pnode):longint;

begin
 If Mode >3 then
  ;
 result:=0;
 if mode=0 then
  exit;
        if (theright^.nodetype=iconstnode) and (theleft^.nodetype=iconstnode) then
         begin
          if mode=3 then
           begin
            result:=2;
             val2:=theright^.value;
             val1:=theleft^.value;
           end
          else
           begin
            result:=1;
            ival2:=theright^.ivalue;
            ival1:=theleft^.Ivalue;
           end;
         end;
        if (theright^.nodetype=constnode) and (theleft^.nodetype=constnode) then
         begin
          result:=2;
          val2:=theright^.value;
          val1:=theleft^.value;
         end;
    if mode>1 then
     begin
       if result=0 then
        begin
         if (theright^.nodetype=constnode) and (theleft^.nodetype=iconstnode) then
         begin
          result:=3;
          val2:=theright^.value;
          val1:=theleft^.ivalue;
         end;
        if (theright^.nodetype=iconstnode) and (theleft^.nodetype=constnode) then
         begin
          result:=4;
          val2:=theright^.ivalue;
          val1:=theleft^.value;
         end;
      end;
    end;
end;

begin
 Result:=0;
 if SimplificationLevel<>0 then
  if param2func then
   result:=DoSimple(SimLevel,expr^.son2left,expr^.son2right)
  else
   result:=DoSimple(SimLevel,expr^.left,expr^.right);

 with expr^ do
  begin
 IF (result>0) and PostProcess then
  begin
   if (operat<>dvdo) then { Divide is special case. If
                                           integer x/y produces a fraction
                                           we want to be able to roll back}
    begin
     if Param2func then
      begin
       dispose(son2right);
       dispose(son2left);
      end
     else
      begin
       dispose(right);
       dispose(left);
      end;
     if result=1 then
      nodetype:=iconstnode
     else
      nodetype:=constnode;
     flags:=[ExprIsConstant];
    end;
   end;
  end;
end;

procedure Checkvarnode(p:pnode);

var treal:arbfloat;
    error:integer;
    tint :Integer;

begin
   TrimLeft(P^.variable);
   TrimRight(p^.variable);
   Val(p^.variable, treal, Error);
   IF (error=0) then {Conversion to real succeeds. Numeric}
    begin
      p^.flags:=[ExprIsConstant];
      if (Pos('.',p^.variable)=0) and (Pos('E',p^.variable)=0) Then
       begin
        Val(p^.variable,tint,Error);
        If error=0 then
         begin
          p^.nodetype:=iconstnode;
          p^.ivalue:=tint;
         end
        else
         begin
          p^.nodetype:=constnode;
          p^.value:=treal;
         end;
       end
      else
       begin
        p^.nodetype:=constnode;
        p^.value:=treal;
       end;
    end;
end;

var tmpval : ArbInt;
    invdummy: pnode;

begin
 case expr^.nodetype of
  VarNode  : CheckVarnode(expr);   {sometimes a numeric value can slip in as constant.
                        (e.g. as people pass it as symbol to taylor or
                        "subst" methods}

  calcnode : begin
              //If not
              internalsimplify(expr^.left); {Reduce left and right as much as possible}
              internalsimplify(expr^.right);
            if isconst(expr^.left) and isconst(expr^.right) then
             begin
              TmpVal:=Setupoperation(expr^.op,SimplificationLevel,true,false);
              if tmpval>0 then
               with expr^ do
               case op of
                addo :
                       if tmpval=1 then
                        ivalue:=ival1+ival2
                       else
                        value:=val1+val2;
               subo  : if tmpval=1 then
                        ivalue:=ival1-ival2
                       else
                        value:=val1-val2;
               mulo  : if tmpval=1 then
                        ivalue:=ival1*ival2
                       else
                        value:=val1*val2;

               dvdo  : begin
                       if tmpval=1 then
                        begin
                         tmpval:=ival1 div ival2;
                         if (tmpval*ival2)=ival1 then {no rounding, OK!}
                          begin
                           Dispose(expr^.right);
                           Dispose(Expr^.left);
                           nodetype:=iconstnode;
                           ivalue:=tmpval;
                          end; {ELSE do nothing}
                        end
                        else
                         begin
                          dispose(expr^.right);
                          dispose(expr^.left);
                          nodetype:=constnode;
                          value:=val1 / val2;
                         end;
                        flags:=[ExprIsConstant];
                        end;
               powo :  If tmpval=1 then
                        begin
                         if ival2<0 then {integer x^-y -> inv (x^y)}
                          begin
                           expr^.nodetype:=funcnode;
                           expr^.son:=NewIConst(IPower(Ival1,-Ival2));
                          end
                         else
                          ivalue:=ipower(ival1,ival2);
                        end
                      else
                        value:=exp(val2*ln(val1));
                 else
                  ExprInternalError(1,ord(Expr^.op));
                end; {case}
            end {if}
           else {At least one node is symbolic, or both types are wrong}
            begin
             With Expr^ do
              if IsConstNil(Left) then
               begin
                Dispose(Left);
               case op of
                addo : begin
                        InvDummy:=Right;
                        Expr^:=Right^;
                        Dispose(InvDummy);
                       end;
                subo: begin
                        invdummy:=right;
                        NodeType:=funcNode;
                        Fun:=Minus;
                        son:=invdummy;
                        Flags:=Son^.Flags;
                      end;
           mulo,powo,dvdo : begin
                        Dispose(Right);
                        nodetype:=IconstNode;
                        ivalue:=0;
                        Flags:=[ExprIsConstant];
                       end;
                   end;
                 end
               else
                if IsConstNil(Right) then
                 begin
                  if expr^.op<>dvdo then {Leave tree for DVD intact because of exception}
                   Dispose(Right);
                 case expr^.op of
       addo,subo : begin
                    InvDummy:=left;
                    Expr^:=left^;
                    Dispose(InvDummy);
                   end;
           mulo  : begin
                    Dispose(Left);
                    nodetype:=IconstNode;
                    Flags:=[ExprIsConstant];
                    ivalue:=0;
                   end;
           powo  : begin
                    Dispose(Left);
                    nodetype:=IconstNode;
                    Flags:=[ExprIsConstant];
                    ivalue:=1;
                   end;
           dvdo  : Raise EDiv0.Create(SExprInvSimp);
             else
                 ExprInternalError(6,ord(Expr^.op));
              end;
            end;
          end;
          With Expr^ Do
           Begin
            IF (nodetype=calcnode) and (Op in [Mulo,Addo]) then
             begin  {Commutative operator rearrangements, move constants to left}
              if (ExprIsConstant IN Right^.flags) and NOT (ExprIsConstant IN Left^.flags) then
              begin
               InvDummy:=Right;
               Right:=Left;
               Left:=InvDummy;
              end;
             IF (right^.nodetype=calcnode) and (right^.Op in [Mulo,Addo]) then
              begin
              end;
             end;
           End;
          end; {case calcnode}

  funcnode: begin
             internalSimplify(expr^.son);
             Case Expr^.fun of
              Minus : if IsConst(expr^.son) then
                       begin
                        InvDummy:=Expr^.Son;
                        expr^:=InvDummy^;
                        if InvDummy^.Nodetype=IconstNode then
                         expr^.ivalue:=-expr^.ivalue
                        else
                         expr^.value:=-expr^.value;
                        dispose(InvDummy);
                       end;
             invx   : begin
                       InvDummy:=Expr^.son;
                       If InvDummy^.nodeType=ConstNode Then
                        begin
                         if InvDummy^.Value=0.0 then
                          Raise EDiv0.Create(SExprInvMsg);
                         Expr^.NodeType:=ConstNode;
                         Expr^.Value:=1/InvDummy^.Value;
                         Dispose(InvDummy);
                        end
                       else
                        if InvDummy^.nodetype=iconstnode then
                         begin
                          if InvDummy^.iValue=0 then
                           Raise EDiv0.Create(SExprinvmsg);
                          If (InvDummy^.iValue=1) or (InvDummy^.iValue=-1) then
                           begin
                            expr^.NodeType:=Iconstnode;
                            Expr^.iValue:=InvDummy^.iValue;
                            Dispose(InvDummy);
                           end;
                         end;
                      end;
                 else {IE check in EvaluateFunction}
                  if (expr^.son^.nodetype=constnode) and (Expr^.fun<>faculx) then {Other functions, only func(real) is simplified}
                   begin
                    val1:=EvaluateFunction(expr^.fun,Expr^.son^.value);
                    dispose(expr^.son);
                    expr^.nodetype:=constnode;
                    expr^.value:=val1;
                   end;
                end; {Case 2}
           end;

  Func2Node : begin
               internalSimplify(expr^.son2left);
               internalSimplify(expr^.son2right);
               case expr^.fun2 of
                powerx : begin
                        TmpVal:=Setupoperation(powo,SimplificationLevel,true,true);
                        if TmpVal>1 then
                         begin
                          If tmpval=1 then
                           begin
                            if ival2<0 then {integer x^-y -> inv (x^y)}
                             begin
                              new(invdummy);
                              invdummy^.nodetype:=iconstnode;
                              invdummy^.ivalue:=ipower(ival1,-ival2);
                              expr^.nodetype:=funcnode;
                              expr^.son:=invdummy;
                             end
                           else
                            expr^.ivalue:=ipower(ival1,ival2);
                           end;
                          end;
                       end;
               stepx : begin
                       {N/I yet}
                       end;
             arctan2x : begin
                         TmpVal:=Setupoperation(powo,SimplificationLevel,false,true);
                         if tmpval>1 then {1 is integer, which we don't do}
                          begin
                           dispose(expr^.right);
                           dispose(expr^.left);
                           expr^.nodetype:=constnode;
                           expr^.value:=arctan2(ival2,ival1);
                          end;
                        end;
             hypotx   :begin
                         TmpVal:=Setupoperation(powo,SimplificationLevel,false,true);
                         if tmpval>1 then {1 is integer, which we don't do}
                          begin
                           dispose(expr^.right);
                           dispose(expr^.left);
                           expr^.nodetype:=constnode;
                           expr^.value:=hypot(ival2,ival1);
                          end;
                        end;
           lognx:      begin
                         TmpVal:=Setupoperation(powo,SimplificationLevel,false,true);
                         if tmpval>1 then {1 is integer, which we don't do}
                          begin
                           dispose(expr^.right);
                           dispose(expr^.left);
                           expr^.nodetype:=constnode;
                           expr^.value:=hypot(ival2,ival1);
                          end;
                        end;
               else
                ExprInternalError(3,ORD(expr^.fun2));
             end;
            end;
{         else
           ExprInternalError(4,ORD(expr^.nodetype));}
       end; {Case 1}
end;

begin
 internalsimplify(exprtree);
 InfixClean:=False; {Maybe changed}
end;

procedure TExpression.SymbolSubst(ToSubst,SubstWith:String);

procedure InternalSubst(expr:Pnode);

begin
 if Expr<>NIL THEN
  case Expr^.NodeType of
   VarNode:   if Expr^.Variable=ToSubst then
               Expr^.Variable:=SubstWith;
   calcnode:  begin
               InternalSubst(Expr^.left);
               InternalSubst(Expr^.right);
              end;
   funcnode:  InternalSubst(Expr^.son);
   func2node: begin
               InternalSubst(Expr^.son2left);
               InternalSubst(Expr^.son2right);
              end;
            end;
end;

begin
 InternalSubst(ExprTree);
end;

function TExpression.SymbolicValueNames:TStringList;

var TheList:TStringList;

procedure InternalSearch(expr:Pnode);

begin
 if Expr<>NIL THEN                      {NIL shouldn't be allowed, and signals corruption. IE? Let it die?}
  case Expr^.NodeType of
   VarNode:  begin
              Expr^.Variable:=UpCase(Expr^.Variable);
              TheList.Add(Expr^.Variable);
             end;
   calcnode:  begin
               InternalSearch(Expr^.left);
               InternalSearch(Expr^.right);
              end;
   funcnode:  InternalSearch(Expr^.son);
   func2node: begin
               InternalSearch(Expr^.son2left);
               InternalSearch(Expr^.son2right);
              end;
            end;
end;

begin
 TheList:=TStringList.Create;
 TheList.Sorted:=TRUE;
 Thelist.Duplicates:=DupIgnore;
 InternalSearch(ExprTree);
 Result:=TheList;
end;

function TExpression.Taylor(Degree:ArbInt;const x,x0:String):TExpression;
{Taylor(x,x0)=sum(m,0,degree, d(f)/d(x))(x0)/ m! * (x-x0)^m)

=   f(x0)+ (x-x0)/1! * df/d(x) + (x-x0)^2  /  2! * d^2(f)/d^2(x) +
       (x-x0)^3  /  3! * d^3(f)/d^3(x) + ....
}

Var TaylorPol    : TExpression;   { The result expression}
    Root,
    Tmp,Tmp2,
    tmp3,tmp4,tmp5 : pnode;       { Always have a nice storage of pointers.
                                    Used to hold all intermediate results}
    I,facul        : ArbInt;      { Loop counters and faculty term}

begin
 TaylorPol:=TExpression.EmptyCreate;      {New expression}
 TaylorPol.ExprTree:=CopyTree(ExprTree);  {make a copy of the parsetree}
 TaylorPol.SymbolSubst(X,X0);             {subst x by x0. All occurances
                                          of  f() refer to x0, not x}
 if Degree>0 then                         {First term only? Or nonsense (negative?)}
                                          {Then ready. 0th term only.}
  begin                                   {Start preparing easy creation of higher terms}
   tmp2:=newcalc(subo,newvar(x),
                      newvar(x0));        {tmp2= x-x0 needed separately for first term}
   tmp4:=Newiconst(5);                    {exponent for x-x0, "a" need to keep a reference}
   tmp3:=newcalc(powo,tmp2,tmp4);         {tmp3= (x-x0)^a}
   tmp5:=newiconst(5);                    {faculty value, "b"=m! also keep a reference for later modification }
   tmp3:=Newcalc(dvdo,tmp3,tmp5);         {tmp3=  (x-x0)^a / b    a and b can be changed}
   facul:=1;                              {Calculate faculty as we go along. Start with 1!=1}
   root:=TaylorPol.ExprTree;              {0th term}
   tmp:=root;                             {term that needs to be differentiated per iteration}
   for i:=1 to Degree do
    begin
     facul:=Facul*i;                      {Next faculty value 1!*1 =1 btw :_)}
     tmp:=intderive(x0,tmp);              {Differentiate f^n(x0) to f^(n+1)(x0)}
     If I=1 then                          {first term is special case. No power }
       tmp2:=NewCalc(mulo,CopyTree(tmp2),tmp) {or faculty needed (^1 and 1!)}
      else
       begin
        tmp5^.Ivalue:=facul;              {Higher terms. Set faculty}
        tmp4^.ivalue:=i;                  {Set exponent (=a) of (x-x0)^a}
        tmp2:=NewCalc(mulo,CopyTree(tmp3),tmp); {multiplicate derivative with (x-x0)^a/b term}
       end;
      root:=NewCalc(addo,root,tmp2);      {String all terms together}
     end;
   DisposeExpr(tmp3);                     {Is only CopyTree()'d, not in new expression}
   TaylorPol.Exprtree:=root;              {Set result}
  end;
 Result:=TaylorPol;
end;

function TExpression.Newton(x:String):TExpression;
{
             f(x)
Newton(x)=x- ----
             f'(x)
}

Var NewtonExpr     : TExpression;   { The result expression}
    Root,
    Tmp,Tmp2,
    tmp3,tmp4,tmp5 : pnode;       { Always have a nice storage of pointers.
                                    Used to hold all intermediate results}
    I,facul        : ArbInt;      { Loop counters and faculty term}

begin
 NewtonExpr:=TExpression.EmptyCreate;      {New expression}

 {Should I test for constant expr here?}

 Tmp:=CopyTree(ExprTree);                 {make a copy of the parsetree
                                            for the f(x) term}
 Tmp2:=intDerive(x,tmp);                  { calc f'(x)}
 Tmp3:=NewVar(x);                         { Create (x)}
 Tmp:=Newcalc(dvdo,tmp,tmp2);             { f(x)/f'(x)}
 tmp:=Newcalc(subo,tmp3,tmp);             { x- f(x)/f'(x)}

 {We built the above expression using a copy of the tree.
     So no pointers into self.exprtree exist. We can now safely assign it
     to the new exprtree}
 NewtonExpr.ExprTree:=tmp;
 NewtonExpr.SimplifyConstants;           {Simplify if f'(x)=constant, and
                                          kill 0*y(x) terms}
 Result:=NewtonExpr;
end;


{
  $Log$
  Revision 1.1  2002/12/15 21:01:26  marco
  Initial revision

}