The Museum of HP Calculators

HP Forum Archive 19

[ Return to Index | Top of Index ]

Symbolic differentiation on HP35s
Message #1 Posted by Peter Niessen on 22 Feb 2009, 9:42 a.m.

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;

}


[ Return to Index | Top of Index ]

Go back to the main exhibit hall