Symbolic differentiation on HP35s

[Peter Niessen]

Hi guys,

I don't know if Vincze reads here anymore. Coming back on how to do symbolic differentiation on the 35s (see this thread), here's a little C++ program which shows how to do symbolic diffrentiation on a RPN expression. You could use the HP35s's storage array as the equivalent of the deque in C++, I guess. Of course, the whole thing would be more academic, since the 16 (?) levels of subroutines won't allow too much of the necessary recursion.

Hoping it's ok to cut and paste here,

// deque.cc

// How to do symbolic differentiation on RPN expressions.

// PN Sat Feb 21 13:15:19 CET 2009

#include <deque>

#include <string>

#include <iostream>

// this is our stack.

// the push_front is equivalent to pressing ENTER

typedef std::deque<std::string> ds_type;

void print (ds_type &d);

ds_type pop_expr (ds_type &d);

ds_type diff_expr (ds_type &d);

ds_type smpl_expr (ds_type &d);

/**********************************************************************/

int main () {

  ds_type expression;

  ds_type diff;

  ds_type smpl;

  expression.push_front ("2");

  expression.push_front ("x");

  expression.push_front ("/");

  std::cout << "This is the orig expression" << std::endl;

  print (expression);

  diff = diff_expr (expression);

  std::cout << "This is the diff'd expression" << std::endl;

  print (diff);

  smpl = smpl_expr (diff);

  std::cout << "This is the smpl'd expression" << std::endl;

  print (smpl);

}

/**********************************************************************/

void print (ds_type &d) {

  std::cout << "***************************************************"

	    << std::endl;

  for (ds_type::const_reverse_iterator i = d.rbegin (); i != d.rend (); i++)

    std::cout << *i << std::endl;

  std::cout << "==================================================="

	    << std::endl;

}

/**********************************************************************/

// returns the next expression from the stack. Modifies the original

// stack! It does so by counting how many operators and operands 

// there

// are. Once it finds that all operators have sufficient operands,

// it

// knows the expression is finished

ds_type pop_expr (ds_type &d) {

  int depth_count = 0;

  int operator_count = 0;

  ds_type result;

  do {

    std::string front = d.front ();

    result.push_back (front);

    d.pop_front ();

    if ("+" == front

	|| "-" ==  front

	|| "*" == front

	|| "/" == front) { // 

			   // binary op, increase by 2, and every

			   // level of operator gets an additional

			   // operator count of 1

      depth_count += 2;

      operator_count += 1;

      continue;

    }

    depth_count -= 1; // it's an operand, decrease by 1

  } while (depth_count - operator_count + 1 > 0);

  return result;

}

/**********************************************************************/

// this is the actual differentiation. It knows the following:

// (f+g)' = f' + g'

// (f-g)' = f' - g'

// (f*g)' = f'*g + f*g'

// (f/g)' = (f'*g-f*g')/g^2

// x' = 1

// (all other)' = 0

ds_type diff_expr (ds_type &d) {

  ds_type result;

  std::string front = d.front ();

  d.pop_front ();

  if ("+" == front || "-" == front) {

    ds_type expr1 = pop_expr (d);

    ds_type expr2 = pop_expr (d);

    ds_type diff1 = diff_expr (expr1);

    ds_type diff2 = diff_expr (expr2);

    // f' +/- g' -> result

    while (diff2.size () > 0) {

      result.push_front (diff2.back ());

      diff2.pop_back ();

    }

    while (diff1.size () > 0) {

      result.push_front (diff1.back ());

      diff1.pop_back ();

    }

    result.push_front (front);

    return result;

  } // f +/- g

  // (f*g)' = f' * g + f * g'

  if ("*" == front) {

    ds_type expr1 = pop_expr (d);

    ds_type expr2 = pop_expr (d);

    // need to make copies because diff_expr removes the thing from

    // the stack

    ds_type undiff1 = expr1;

    ds_type undiff2 = expr2;

    ds_type diff1 = diff_expr (expr1);

    ds_type diff2 = diff_expr (expr2);

    // f' * g -> result

    while (diff1.size () > 0) {

      result.push_front (diff1.back ());

      diff1.pop_back ();

    }

    while (undiff2.size () > 0) {

      result.push_front (undiff2.back ());

      undiff2.pop_back ();

    }

    result.push_front ("*");

    // f * g' -> result

    while (undiff1.size () > 0) {

      result.push_front (undiff1.back ());

      undiff1.pop_back ();

    }

    while (diff2.size () > 0) {

      result.push_front (diff2.back ());

      diff2.pop_back ();

    }

    result.push_front ("*");

    // f' * g + f * g'

    result.push_front ("+");

    return result;

  } // f * g

  // (f/g)' = (f'g - fg')/g^2

  if ("/" == front) {

    ds_type expr1 = pop_expr (d);

    ds_type expr2 = pop_expr (d);

    ds_type undiff1 = expr1;

    ds_type undiff11 = expr1;

    ds_type undiff12 = expr1;

    ds_type undiff2 = expr2;

    ds_type diff1 = diff_expr (expr1);

    ds_type diff2 = diff_expr (expr2);

    // f * g'

    while (diff2.size () > 0) {

      result.push_front (diff2.back ());

      diff2.pop_back ();

    }

    while (undiff1.size () > 0) {

      result.push_front (undiff1.back ());

      undiff1.pop_back ();

    }

    result.push_front ("*");

    // f' * g

    while (undiff2.size () > 0) {

      result.push_front (undiff2.back ());

      undiff2.pop_back ();

    }

    while (diff1.size () > 0) {

      result.push_front (diff1.back ());

      diff1.pop_back ();

    }

    result.push_front ("*");

    // f' * g - f * g'

    result.push_front ("-");

    

    while (undiff11.size () > 0) {

      result.push_front (undiff11.back ());

      undiff11.pop_back ();

    }

    while (undiff12.size () > 0) {

      result.push_front (undiff12.back ());

      undiff12.pop_back ();

    }

    result.push_front ("*");

    result.push_front ("/");

    return result;

  } // f/g



  // this one's easy

  if ("x" == front) {

    result.push_front ("1");

    return result;

  }

  // if we got here, it must be a constant

  result.push_front ("0");

  return result;

    

}

/**********************************************************************/

// simplify the expression by using rules:

// x + 0 = x

// x - 0 = x

// x * 1 = x

// x * 0 = 0

// x / 1 = x

// 0 / x = 0

ds_type smpl_expr (ds_type &d) {

  ds_type original = d;

  // check what type of expression it is

  std::string front = d.front ();

  d.pop_front ();

  if (front == "+") {

    ds_type expr1 = pop_expr (d);

    ds_type smpl1 = smpl_expr (expr1);

    ds_type expr2 = pop_expr (d);

    ds_type smpl2 = smpl_expr (expr2);

    // now simplify: if expr1 is 0, just return expr2

    if (smpl1.size () == 1 && smpl1.front () == "0")

      return smpl2;

    if (smpl2.size () == 1 && smpl2.front () == "0")

      return smpl1;

    // if the rules above don't apply, return the sum of the two

    // simplified operands

    ds_type result;

    while (smpl1.size () > 0) {

      result.push_front (smpl1.back ());

      smpl1.pop_back ();

    }

    while (smpl2.size () > 0) {

      result.push_front (smpl2.back ());

      smpl2.pop_back ();

    }

    result.push_front ("+");

    return result;

    

  }

  if (front == "-") {

    ds_type expr1 = pop_expr (d);

    ds_type smpl1 = smpl_expr (expr1);

    ds_type expr2 = pop_expr (d);

    ds_type smpl2 = smpl_expr (expr2);

    if (smpl1.size () == 1 && smpl1.front () == "0")

      return smpl2;

    // if the rules above don't apply, return the difference of the two

    // simplified operands

    ds_type result;

    while (smpl2.size () > 0) {

      result.push_front (smpl2.back ());

      smpl2.pop_back ();

    }

    while (smpl1.size () > 0) {

      result.push_front (smpl1.back ());

      smpl1.pop_back ();

    }

    result.push_front ("-");

    return result;

  }

  if (front == "*") {

    ds_type expr1 = pop_expr (d);

    ds_type smpl1 = smpl_expr (expr1);

    ds_type expr2 = pop_expr (d);

    ds_type smpl2 = smpl_expr (expr2);

    // now simplify: if expr1 is 1, just return expr2

    if (smpl1.size () == 1 && smpl1.front () == "1")

      return smpl2;

    if (smpl2.size () == 1 && smpl2.front () == "1")

      return smpl1;

    // x * 0 = 0

    if (smpl1.size () == 1 && smpl1.front () == "0"

	|| smpl2.size () == 1 && smpl2.front () == "0") {

      ds_type result;

      result.push_front ("0");

      return result;

    }

    // if the rules above don't apply, return the product of the two

    // simplified operands

    ds_type result;

    while (smpl1.size () > 0) {

      result.push_front (smpl1.back ());

      smpl1.pop_back ();

    }

    while (smpl2.size () > 0) {

      result.push_front (smpl2.back ());

      smpl2.pop_back ();

    }

    result.push_front ("*");

    return result;

  }

  if (front == "/") {

    ds_type expr1 = pop_expr (d);

    ds_type smpl1 = smpl_expr (expr1);

    ds_type expr2 = pop_expr (d);

    ds_type smpl2 = smpl_expr (expr2);

    if (smpl1.size () == 1 && smpl1.front () == "1")

      return smpl2;

    if (smpl2.size () == 1 && smpl2.front () == "0") {

      ds_type result;

      result.push_front ("0");

      return result;

    }

    // if the rules above don't apply, return the ratio of the two

    // simplified operands

    ds_type result;

    while (smpl2.size () > 0) {

      result.push_front (smpl2.back ());

      smpl2.pop_back ();

    }

    while (smpl1.size () > 0) {

      result.push_front (smpl1.back ());

      smpl1.pop_back ();

    }

    result.push_front ("/");

    return result;

  }

  // we got here, there was nothing so simplify. So just return the

  // original expression

  return original;

}