Symbolic differentiation on HP35s



#2

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;

}


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP35s Program Four Slings Lift Calculation Jean-Marc Biram (Australia) 2 1,209 12-16-2013, 07:21 PM
Last Post: Jean-Marc Biram (Australia)
  HP35s Calculator Max Rope Tension Program Jean-Marc Biram (Australia) 10 2,450 12-12-2013, 12:03 AM
Last Post: Jean-Marc Biram (Australia)
  HP Prime suggestion to avoid Numeric/Symbolic confusion Chris Pem10 4 1,090 11-19-2013, 05:49 AM
Last Post: bluesun08
  [HP-Prime] AMBIGUITY between Numerical Calculation (HOME) and Numerical/Symbolic Calculation (CAS mode) CompSystems 2 833 08-18-2013, 07:06 PM
Last Post: CompSystems
  Trouble entering a HP35s program line Arno 2 852 04-05-2013, 06:28 PM
Last Post: Arno
  HP35s scientific calculator GREG W THOMAS 4 1,073 03-22-2013, 06:49 AM
Last Post: Thomas Radtke
  HP35s "MEMORY CLEAR" flashes Mark Paris 1 723 08-31-2012, 07:35 PM
Last Post: Bart (UK)
  Symbolic limit of a function of 2 variables, HP49/50 Gilles Carpentier 0 509 08-26-2012, 10:28 AM
Last Post: Gilles Carpentier
  HP35S keyboard Nick R 8 1,577 08-01-2012, 01:27 PM
Last Post: Dave Shaffer (Arizona)
  Where should I post new HP35s bugs Andres Capdevila 33 4,650 03-13-2012, 01:00 PM
Last Post: Jeff O.

Forum Jump: