Symbolic derivation


symbolic derivation
deadline online
how to find roots
numerical methods
Newton method
bisection method
plot graphs
support center

Problem from ACM Southeastern European Regional Programming Contest
(Bucharest, Romania - October 21, 2000).

Input File: standard input

Program Source File: derivation.cpp

Write a program that performs symbolic derivation f'(x) = df(x)/dx of a given function f(x). The function f(x) is defined by an expression which may contain the following operations: + (addition), - (subtraction), * (multiplication), / (division), and ln (natural logarithm). Besides, the operands may be the variable x and numerical constants. The expression may contain arbitrarily nested sub-expressions in parentheses ( ). The expression is given in a usual, infix form, such as:

(2*ln(x+1.7)-x*x)/((-7)+3.2*x*x)+(x+3*x)*x

Numerical constants have the form d.d, with an optional sign (+ or -), where the number of digits both in integer and decimal parts are arbitrary. The input expression is guaranteed to be correct (no syntax error can occur).

The output expression should be written in infix form. It should not be optimized for human reading. This means, it can contain redundancies, such as 0*x, 1*x, 0+x, etc. The derivation should be performed using the following rules:

1. The operators * and / are of higher priority than the operators + and -. Parentheses may change the priorities as usually.

2. The operators +, -, *, and / are left-associative, meaning that they group from left to right: a*b*c = (a*b)*c, a/b/c = (a/b)/c, a/b*c = (a/b)*c, etc.

3. The rules for derivation are:

(a + b)' = a' + b'

(a - b)' = a' - b'

(a * b)' = (a' * b + a * b')

(a / b)' = (a' * b - a * b') / b^2 Note: use b^2 and not (b*b) for presentation

ln(a)' = (a')/(a)

x' = 1

const' = 0

4. While producing the symbolic derivation, use parentheses for output strictly as stated in the previous rule. Do not perform presentation optimizations, such as 0*a = 0, 1*a = a, etc.

The input is a textual file which has one f(x) definition per line. The input lines do not have blanks. The output should contain lines with corresponding symbolic derivations f'=df/dx, one line for each f. The strings representing f(x) and f'(x) are guaranteed to have no more than 100 characters.

code.h


// Copyright 2004 Ionut Alex. Chitu
// http://deadline.3x.ro

#ifndef CODE_H 
#define CODE_H 

#include <string> 

using namespace std;

enum Coduri { 
     CodNumar=0, CodPd=1, CodPi=2, CodOperator=30, 
     CodPlus=30, CodMinus=31, CodInmultit=32,
     CodImpartit=33, CodMinusUnar=34, CodPutere=35, 
     CodFunctie=101, CodLn=101, CodX=102, CodNec=120 };

enum Tipuri { 
     TipNumar, TipPd, TipPi, TipConst, TipFunctie, 
     TipOpBinar, TipOpUnar, TipNec };

typedef long double ldouble;
typedef unsigned char byte;

typedef struct element {
  byte info, tip, para;
  string num;
  element (byte a=CodNec, byte b=TipNec) { info=a; tip=b; para=0; }
  element (string a) { info=CodNumar; tip=TipNumar; num=a; para=0; }
} element;

const element ElemNec(CodNec, TipNec);
const element ElemPlus(CodPlus, TipOpBinar);
const element ElemMinus(CodMinus, TipOpBinar);
const element ElemInmultit(CodInmultit, TipOpBinar);
const element ElemImpartit(CodImpartit, TipOpBinar);
const element ElemPutere(CodPutere, TipOpBinar);
const element ElemMinusUnar(CodMinusUnar, TipOpUnar);
const element ElemPd(CodPd, TipPd);
const element ElemPi(CodPi, TipPi);
const element ElemLn(CodLn, TipFunctie);
const element ElemX(CodX, TipConst);
      
// 0=+, 1=-, 2=*, 3=/, 4=unary minus, 5=^
const byte prio1[6][6] =
            { {0, 0, 0, 0, 0, 0}, //+
              {0, 0, 0, 0, 0, 0}, //-
              {1, 1, 0, 0, 0, 0}, //*
              {1, 1, 0, 0, 0, 0}, // /
              {0, 0, 0, 0, 0, 0}, //- unary
              {1, 1, 1, 1, 1, 0} }; // ^

const byte prio2[6][6] =
            { {0, 0, 0, 0, 0, 0}, //+
              {1, 1, 0, 0, 1, 0}, //-
              {1, 1, 0, 0, 1, 0}, //*
              {1, 1, 0, 1, 1, 0}, // /
              {1, 1, 0, 0, 1, 0}, //- unary
              {1, 1, 1, 1, 1, 0} }; // ^
              
const byte ord[6][6] =
            { {1, 1, 0, 0, 0, 0}, //+
              {1, 1, 0, 0, 0, 0}, //-
              {1, 1, 1, 1, 0, 0}, //*
              {1, 1, 1, 1, 0, 0}, // /
              {1, 1, 0, 0, 1, 0}, //- unary
              {1, 1, 1, 1, 1, 1} }; // ^
                            		      		      
element decod (const string& s)
{
    if (s.compare ("ln")) return ElemX;
    else return ElemLn;
}

void GetNext (const string& s, element &last, int &i) // get next token
{
   char c=s[i];
   if ( (c>='0' && c<='9') || c=='.') {
       int k=i;
       for (i++; (s[i]>='0' && s[i]<='9') || s[i]=='.'; i++);
       last=element(s.substr(k, i-k)); i--;
     }
   else if (c>='a' && c<='z') {
       int k=i;  
       for (i++; s[i]>='a' && s[i]<='z'; i++);
       last=decod(s.substr(k, i-k)); i--;
   }
   else
   switch (c) {
       case '+': last=ElemPlus; break;
       case '-': 
       if (last.tip==TipPd || last.tip==TipOpBinar) {
    	  if ((s[i+1]>='0' && s[i+1]<='9') || (s[i+1]=='.')) {
    	      int k=i;
              for (i++; (s[i]>='0' && s[i]<='9') || s[i]=='.' ; i++);
    	      last=element(s.substr(k, i-k));
    	    }		  
	      else last=ElemMinusUnar;
	      i--;
	    }    
       else last=ElemMinus;
       break;
       case '/': last=ElemImpartit; break;
       case '*': last=ElemInmultit; break;
       case '^': last=ElemPutere; break;
       case '(': last=ElemPd; break;
       case ')': last=ElemPi; break;
   } //switch
   i++;

}

#endif

derivation.cpp


// http://deadline.3x.ro

#include <stack> 
#include <iostream> 
#include "code.h" 

using namespace std;

typedef struct Tree 
{
  element lexem;
  Tree *left, *right;
  Tree () { lexem=ElemNec; }
  Tree (element el):lexem(el) { }     
} tree;

typedef tree* ptree;

class symbolic
{
   public :
   symbolic(string s);   
   inline string derive_tree();   
   private: 
   ptree t;
   string display (ptree t);
   string derive(ptree t);
};

inline void add2tree (stack<ptree> &cap, element &crt, byte p=0)
{
   crt.para=p;
   ptree tmp = new tree(crt);
   if (crt.tip==TipFunctie || crt.tip==TipOpUnar) {
        tmp->left = cap.top(); cap.pop();
      }    	    
   else if (crt.tip==TipOpBinar) {
        tmp->right = cap.top(); cap.pop();
        tmp->left = cap.top(); cap.pop();
	  }
   cap.push (tmp);
}

symbolic::symbolic (string s)
{
   s="("+s+")";     
   stack<element> top;
   stack<ptree> cap;
   element cod, cod1, crt=ElemNec;
   for(int i=0, l=s.length(); i<l;) {
       GetNext (s, crt, i);
       if (crt.tip==TipPd || crt.tip==TipFunctie || crt.tip==TipOpUnar)
    	 top.push (crt);	 
       else if (crt.tip==TipPi) {
        	 cod=top.top(); top.pop();
        	 if (cod.tip==TipPd) { 
                 if (top.empty()) continue;       	     
        	     cod1=top.top();
         	     if (i<l)
                   if (cod1.tip!=TipFunctie)  { // tree stack
                     if (cap.empty()) continue;
                     ptree tmp=cap.top();                       
        	         cap.pop();
        	         tmp->lexem.para=1;
        	         cap.push(tmp);
        	       }    
        	       else {
        	         top.pop();
        	         add2tree (cap, cod1);  
        	       }    
        	 }          
        	 else while (cod.tip!=TipPd) {
                cod1=top.top(); top.pop();  
        	    if (cod1.tip!=TipPd) add2tree (cap, cod);
        	    else { 
                 if (!top.empty()) {                         
        	        cod1=top.top();
                    if (cod1.tip==TipFunctie) {
                      top.pop();                        
                      add2tree (cap, cod);                  
            	      add2tree (cap, cod1);
            	     }    
            	    else
                      add2tree (cap, cod, (byte)(i<l)); 
                  }
                  else
                      add2tree (cap, cod, (byte)(i<l));      
                    break;  
            	  }        	         
        	    cod=cod1;
        	 }    	 
            }
       else if (crt.tip==TipOpBinar) {
    	   while (!top.empty()) {
                cod=top.top();
                if ((cod.tip==TipOpBinar || cod.tip==TipOpUnar) && 
                     ord[cod.info-CodOperator][crt.info-CodOperator]) {
                     top.pop();
       		         add2tree (cap, cod);
                  }
                else break;
    		  }
    	   top.push (crt);
	   }
	   else add2tree (cap, crt);
     }
   // create tree
   t=cap.top(); cap.pop();
}

string symbolic::derive (ptree t)
{
  string s, v;
  if (t) {
      if (t->lexem.para) s+="(";
      switch (t->lexem.info) {
        case CodX: 
            s+="1"; 
            break;
        case CodNumar: 
            s+="0"; 
            break;
        case CodPlus:
            s+=derive(t->left);        
            v=derive(t->right);
            if (v[0]!='-') v="+"+v;
            s+=v;
            break;
        case CodMinus: 
            s+=derive(t->left);
            v=derive(t->right);
            if (v[0]=='-') v[0]='+';
            else if (v[0]=='+') v[0]='-';
            else v="-"+v;
            s+=v;                    
            break;
        case CodMinusUnar: 
            v=derive(t->left);
            if (v[0]=='-') v[0]='+';
            else v="-"+v;
            s+=v;                         
            break;
        case CodInmultit: 
            s+="(";
            s+=derive(t->left);        
            s+="*";
            s+=display(t->right);
            v=display(t->left);
            if (v[0]!='-') v="+"+v;
            s+=v;
            s+="*";
            s+=derive(t->right);        
            s+=")";
            break;
        case CodImpartit: 
            s+="(";
            s+=derive(t->left);        
            s+="*";
            s+=display(t->right);
            v=display(t->left);
            if (v[0]=='-') v[0]='+';
            else v="-"+v;
            s+=v;            
            s+="*";
            s+=derive(t->right);                
            s+=")/";
            s+=display(t->right);
            s+="^2";        
            break;
        case CodLn:
            s+="(";
            s+=derive(t->left);        
            s+=")/(";
            s+=display(t->left);
            s+=")";
     }   
      if (t->lexem.para) s+=")";
} // if
  return s;
}

string symbolic::display (ptree t)
{
  string buffer;
  if (t) {
      if (t->lexem.para) buffer+="(";
      byte cod=t->lexem.tip;
      if (cod==TipFunctie)
	       buffer+="ln("+display(t->left)+")";
      else if (cod==TipNumar)
           buffer+=t->lexem.num;      
      else if (cod==TipOpBinar) {
        buffer+=display(t->left);
        switch (t->lexem.info) {
          case CodPlus     : buffer+="+"; break;
          case CodMinus    : buffer+="-"; break;
          case CodInmultit : buffer+="*"; break;
          case CodImpartit : buffer+="/"; break;
         }  
        buffer+=display(t->right);
       }
      else if (cod==TipOpUnar) buffer+="-"+display(t->left);
      else buffer+="x";
      if (t->lexem.para) buffer+=")";
   }
  return buffer;
}

inline string symbolic::derive_tree () 
{ 
   return derive(t);
}

int main () 
{
  string s;
  while(cin>>s) {
    symbolic a(s);
    cout<<a.derive_tree()<<endl;
  }    
}

Syntax highlighting by WebCpp.

Evaluate functions easily using Microsoft Script Control.

DeadLine OnLine - free equation solver. Copyright 2003-2007 Ionut Alex. Chitu. | Contact | Sitemap