Implementing Symbolic Algebra & Calculus in Perl

Jonathan is currently working on an undergraduate degree in Computer Science at the University of Cambridge. He can be contacted by e-mail at [email protected].

Computers are good at working with numbers. With the ability to perform millions of calculations in a fraction of a second, complex problems can be solved numerically. Sometimes, however, it's useful to find an analytical solution to a problem, or at least part of a problem, that is initially expressed analytically.

One advantage is that a solution is obtained without any loss of accuracy. A solution computed from a numerical method may suffer from this, even if the method is mathematically supposed to give an exact answer due to the limits of floating-point arithmetic, for example.

Another advantage is that the more general solution provided by an analytical method may make further computation cheaper. Suppose, for example, that an analytical solution is found for the derivative of an expression. If the derivative needs to be evaluated at many points, evaluating the derivative at each point simply means substituting values for the appropriate variable.

This article will look at a way of representing mathematical expressions so they can easily be manipulated and evaluated. In part, it will reference three modules I have worked on recently, which may or may not be a good starting point for working with mathematical expressions analytically in Perl. These modules are, namely:

• Math::Calculus::Expression, a module for representing an expression.
• Math::Calculus::Differentiate, a module that computes the (analytical) derivative of an expression with respect to a single variable.
• Math::Calculus::NewtonRaphson, which demonstrates a conjunction of numerical and analytical methods, and provides a way of solving an equation, currently in just a single variable.

These modules are all available on CPAN, and can be installed in the usual way:

perl -MCPAN -e 'install Math::Calculus::Expression' 
perl -MCPAN -e 'install Math::Calculus::Differentiate'
perl -MCPAN -e 'install Math::Calculus::NewtonRaphson'  

Representing Expressions

An expression in human-readable form, such as sin(2*x), is not easy to manipulate in a computer program. Regular expressions are not suitable because the language of valid mathematical expressions is not regular. Perl 5's regexes are more powerful than real regular expressions and a recursively defined regex could help here, but this probably would not have led to an easy way of implementing programs to manipulate expressions. Therefore, having a separate internal representation seemed to be the way forward.

The requirements for the representation were simple: It needed to be able to represent a mathematical expression, be relatively easy to translate to and from a human-readable form, and be simple to write programs to operate on. After some pondering, the possibility of using a binary tree-based structure came to mind. Under this scheme, a node can contain any of the following:

• An operator, such as +, -, *, / (divide), and ^ (raise to a power). These would always have two branches—one for each operand.
• A function, such as sin, cos, tan, ln, exp, and the like. These would always have one branch, as these functions only take one parameter.
• A variable or constant, such as x, k, 2.4, and the like. These would have no branches.

Any branch in the tree can be assumed to have brackets around it without any loss of expression. This is because an expression such as 5*x+2 is equal to the expression (5*x)+(2), which is also equal to ((5)*(x))+(2). This suggests that this representation is capable of representing all mathematical expressions, for they can all be translated into an "everything bracketed" form.

Out of this analysis comes a simple method of converting an expression in human-readable form to a binary tree. Taking the expression 5*x+2 as an example, it is clear that the tree built from it should look like Figure 1.

Notice that the left subtree is the same tree that would be found if the expression 5*x was to be turned into a tree, and that the right node is the same output that would be expected if the constant 2 was put into the tree builder. This is a strong hint that recursion will be useful here. In fact, the entire process of building the tree can be summarized in three steps. First, find the operator to use for this node. Second, recursively build the left subtree, passing everything to the left of our chosen operator. Third, recursively build the right subtree, passing everything to the right of our chosen operator.

The second and third steps are simple, but the first requires a little thought. It is important that the correct operator is chosen at each step, and for this the tree builder needs to understand precedence. To help keep the tree building function itself neat, another function can be defined that takes two operators and returns True if the precedence of the first operator is lower than that of the second. This means that finding the operator to choose to split on is a case of scanning through the expression character by character looking for operators that are not nested inside any brackets, then when one is found, checking if it has lower precedence than the operator that is currently regarded as the best one to choose. If it does have lower precedence, it is marked as the current best choice operator, and the scan continues to the end of the expression.

This isn't quite the whole story; there are two special cases to be aware of. These should be checked for before the previously described routine is entered, for if either are true, it would fail.

The first of these special cases involves the expression to build a tree from a constant or variable. If this is the case, it will either be a number (which may or may not be an integer) or a single letter. In both cases, a - sign may precede it. This case can be identified using a simple pattern match; /^-?(\d+(.\d|[A-Za-z])$/ will work.

The second special case is where the entire expression is a function of a subexpression. In this case, the expression will start with the name of a recognized function, such as sin, followed by an opening bracket. Again, a - sign may precede it. The name of the function is stored where the operator would otherwise be located, and the left branch of the tree is generated by recursively calling the tree builder on the function's parameter. There is no right branch.

Converting the tree representation back into a human-readable form can also be achieved through a (much simpler) recursive process. A pretty printing function can be defined as follows:

• If the current node is a constant or variable, this can simply be returned.
• If the current node is a function, recursively call the pretty print function on the left branch of the tree, then return something of the form "<function name>(<left branch>)."
• If the current node is an operator, recursively call the pretty print function on the left branch and the right branch, then return something of the form "(<left branch>) <operator> (<right branch>)."

Some extra tracking of operators can be done to work out whether brackets are needed or if they can be omitted because precedence rules make them redundant.

So far, two of the three requirements for the representation are satisfied; it can represent mathematical functions and can easily and fairly cheaply be converted to and from a human-readable representation. The third condition was that it should be easy to manipulate mathematical expressions using this representation. The binary tree-based structure has some useful properties. It is now trivial to identify and extract subexpressions; every subtree is a valid expression. Precedence is completely determined by depth in the tree, and there is no need to do any more work to check operator precedence when manipulating the expression. If this representation has been formed, it also means that the expression is valid, removing much of the responsibility for error checking from code manipulating the expression. The general benefits of separating parsing from evaluation are gained, too: Now there is just one place to fix parsing bugs.

Math::Calculus::Expression

The Math::Calculus::Expression module implements a class that represents a mathematical expression using the binary tree method described earlier. Remember that an expression is not an equation; there must be no "=" symbols or other statements about mathematical equivalence. It understands the operators +, -, *, / (divide) and ^ (raise to power) and their precedence, as well as the functions sin, cos, tan, sec, cosec, cot, their inverses (add an a to the start of those names—for example, asin is inverse sin) and hyperbolic equivalents (add an h to the end—for example, sech is hyperbolic sec), ln (natural logarithm), and exp (base of natural logarithms to the power of the parameter). Use of brackets to control precedence is supported. All subexpressions must be separated by an operator, so 5*x and (x+2)*(x-3) are valid but 5x and (x+2)(x-3) are not. Spacing does not matter, and most certainly does not control precedence. 2+x * 5 means exactly the same as 2 + x*5.

A new expression is instantiated using the new method of the Expression class:

use Math::Calculus::Expression;
my $expr = Math::Calculus::Expression->new();

Various methods can then be called on the expression object. The setExpression method takes an expression in human-readable form as its only parameter and attempts to convert it into the internal representation. If this fails, which will happen if the expression is invalid, then a False value will be returned. A traceback that may contain useful information is also available by calling the getTraceback method, and a short error message can be obtained by calling the getError message:

unless ($expr->setExpression(
         '2*sin(x^2 + 4) + cosh(x)')) {
  print "Error: " . $expr->getError . "\n";
  print "Traceback: \n" . $expr->getTraceback;
}

Expressions stored in an expression object can be retrieved in a human-readable form using the getExpression method:

# Prints 2*sin(x^2 + 4) + cosh(x)
print $expr->getExpression;

Note that getExpression may not return exactly the same string that was provided when setExpression was called. The internal representation retains the meaning of the expression, but not spacing and extraneous brackets.

It is often useful to tell the module which letters correspond to variables; any that don't will be assumed to be constants. The addVariable method is used to declare a letter as a variable.

$expr->addVariable('x');

The simplify method attempts to make the expression simpler. Its bag of tricks include constant evaluation (2*3*x becomes 6*x), redundant-term elimination (0*sin(x) and 2*x-2*x will become 0, and if part of a more complex expression, will disappear), and null operator elimination (x^1 and 1*x and x+0 will all become x). The method is simply invoked on an expression object and modifies the expression that it stores.

$expr->simplify;

The evaluate method takes a hash, mapping any variables and named constants (collectively letters) in the expression to numerical values, and attempts to evaluate the expression and return a numerical value. It fails if it finds letters that have no mapping or if an error such as division by zero occurs. The currently stored expression is not modified, so this method can be called to evaluate the expression at many points.

$value = $expr->evaluate(x => 0.5);

# Use defined as 0 is false but valid
unless (defined($value)) { 
  print "Error: " . $expr->getError . "\n";
}

The sameRepresentation method takes another expression object as its parameter and returns True if that expression has the same internal representation as the expression the method is invoked on; under the current implementation, that means identical trees. Be careful—while it can be said that if two expressions have the same representation, they are equal; it would be false to say that if they have different representations, they are not equal. It is clear that x+2 and 2+x are equal, but their internal representation may well differ. As the internal representation and its semantics may change over time, the only thing that is safe to say about this method is that if it returns True, the expressions are equal. If it returns False, nothing can be said about equality.

print "Same representation" if $expr-> 	sameRepresentation($expr2);

Finally, the clone method returns a deep copy of the expression object ("deep copy" meaning that if the original is modified, the copy will not be affected and vice versa):

$exprCopy = $expr->clone;

Operating on an Expression

The Math::Calculus::Expression class uses a hash for each node. The hash has keys operator (storing the operator or function name), operand1 (for the left branch), and operand2 (for the right branch). The values corresponding to the branch keys (operand1 and operand2) may be undefined, contain a number or letter, or contain a reference to another hash. Using a hash-based structure helps with program readability and allows others modules that need to do more sophisticated manipulation to augment the tree with their own data by simply putting extra entries into the hash.

Generally, tree operations are implemented using recursion, as a depth-first exploration of the tree is usually required. A breadth-first exploration may also be useful in some cases, and this could also be implemented. This article will take a look at the depth-first approach. Operations can be subdivided into three categories: those that modify the tree, those that use this tree to build a new tree, and those that use the tree to generate something else (such as a single value). Knowing when to use the third of these is fairly easy, but it's not always obvious when choosing between the first and the second. A good rule of thumb is to consider whether the operation changes the values or the existing structure of the tree. If the structure is changing, then it will usually be simpler to build a new tree.

As a simple example, consider writing a module that can substitute an expression or a constant for a particular letter. A good starting point is to inherit from Math::Calculus::Expression, as this will provide a constructor and other useful methods:

package Math::Calculus::Substitute;
use strict;
use Math::Calculus::Expression;
our @ISA = qw/Math::Calculus::Expression/;

The structure of the expression tree may change if a complex expression is substituted for a letter, but note that this is simply an extension and not a change to the structure that is already in place—basically, a single value is being changed, but it may be changed to be a subtree rather than another value. Therefore, this problem is simple enough that it's best to modify the existing tree. A procedure can be written that will recurse down the tree looking for a branch that is equal to the letter that is being sought and replace it with whatever replacement has been provided, which it will assume is ready to be substituted. This is OK because it will be a private method rather than the public interface to the module.

sub recSubstitute {
  # Grab the parameters.
  my ($self, $tree, $letter, $sub) = @_;
  
  # Look at each branch.
  my $branch;
  foreach $branch ('operand1', 'operand2') {
    # If it's a reference to a subtree, recurse.
    if (ref $tree->{$branch}) {
      $self->recSubstitute($tree->{$branch}, 
                           $letter, $sub);
    
    # Otherwise, if it's the letter...
    } elsif ($tree->{$branch} eq $letter) {
      # Substitute.
      $tree->{$branch} = $sub;
    }
  }
}

Finally, the public method needs to be implemented. It needs to use a couple of the private methods of the Expression class that exist for the usage of subclasses, namely getExpressionTree to get a copy of the raw expression tree and setError to set or clear the current error message. As the substitution could be a letter/number or an expression object (but not a raw tree—remember that the tree representation is internal), this method needs to ensure that what is passed to recSubstitute is something that can be substituted straight into the tree.

sub substitute {
  # Grab parameters.
  my ($self, $letter, $sub) = @_;

  # Check letter really is a letter.
  if ($letter !~ /^[A-Za-z]$/) {
    $self->setError('A valid letter must 
                             be passed.');
    return undef;
  }
  
  # If the sub is a reference, then it 
  # should be pointing to an Expression
  # object. Extract the raw expression.
  my $realSub;
  if (ref($sub)) {
    $realSub = $sub->getExpressionTree;
  } else {
    $realSub = $sub;
  }
  
  # Get a reference to this objects's 
  # expression tree.
  my $tree = $self->getExpressionTree;
  
  # Make the recursive call.
  $self->recSubstitute($tree, $letter, $realSub);
  
  # Clear the error and return success.
  $self->setError('');
  return 1;
}

The next two sections will take a quick look at a couple of other modules that inherit from the Expression class and do some more complex manipulations.

Math::Calculus::Differentiate

Differentiation is one part of the calculus. It involves taking a function and finding its gradient function—another function that defines how steep the function is at any given point. In many cases, the function to differentiate can be written in the form y=f(x). As f(x) is an expression and finding the derivative is a clearly defined set of manipulations on the expression, the tools described so far provide all that is needed to implement a module to do differentiation analytically.

As many of the manipulations involve significant structural changes to the existing tree, the Differentiate module builds a new tree as it walks the existing one. This is hidden from the end user of the module as the public differentiate method takes the new tree and sets the current expression object to point to it. This is a good thing to do in general, as it leaves the user with only two kinds of operations to worry about: those that don't change the tree and return a result, and those that do modify the tree and return success or failure.

Usage of the module is simple: Call the differentiate method and pass in the variable to differentiate with respect to. Note that all differentiation done by this module is partial; if there are multiple variables, then others are taken to be constants, as is the norm in partial differentiation. The following example shows the module at work:

use Math::Calculus::Differentiate;

# Create new expression object.
my $expr = Math::Calculus::Differentiate->new;

# Set expression and variable.
$expr->setExpression('x^2 + sin(x)');
$expr->addVariable('x');

# Differentiate with respect to x and simplify.
$expr->differentiate('x');
$expr->simplify;

# Output result.
print $expr->getExpression; # Prints 2*x + cos(x)

Note that error checking has been omitted from this example to keep it simple, and should of course be included in any real program that uses this module.

Math::Calculus::NewtonRaphson

Solving equations analytically is sometimes impossible for people, let alone computers. This is a case where a numerical method is almost always the right thing to use. The Newton Raphson method of solving an equation requires that it is rearranged into the form f(x)=0. Again, f(x) is an expression and it is therefore suitable to use the system described in this article to manipulate it.

The Newton Raphson method is iterative; that is, it takes an initial guess and produces a solution. That solution can be fed back in to produce a better solution. This process can be repeated until a solution has been found to the desired accuracy. Being able to use the Newton Raphson method involves knowing the derivative of f(x), f(x) itself, and having a good initial guess. The first of these can be calculated with the Math::Calculus::Differentiate module, and the second two are things that the user of the module supplies.

This is an example of a module that produces a result without modifying the existing expression—at least, it does from the user's point of view. Under the hood, it will need to make a copy of the expression and differentiate it, then form a new expression involving the original one and the derivative. This can then be passed to the evaluate method inherited from the Expression class. A simple example of the module at work follows:

use Math::Calculus::NewtonRaphson;

# Create new expression object.
my $expr = Math::Calculus::NewtonRaphson;

# Set expression and variable.
$expr->setExpression('x^2 - 4');
$expr->addVariable('x');

# Solve with Newton Raphson - initial guess of 3.
my $result = $expr->newtonRaphson(3);

# Display result, which will 
# be 2 as 2^2 - 4 = 4 - 4 = 0
print $result; # Prints 2

It is possible that a solution will not be found; in which case, the newtonRaphson method will return undef. Any real-world program should check for this case.

Where From Here

There is still plenty of scope for improving the simplifier and writing more modules. One big problem with what I have described in this article is that there is currently no way to test for equality—whether two expressions have the same representation in the current tree format is too tight a condition. One solution is to try to find a canonical representation—one in which two things that are equal will always get the same representation. This is what happens in binary decision diagrams, and it enables us to cheaply determine equality by comparing representations. One problem with this is that work may need to be done to ensure that the expression stays in its canonical representation, and enabling easy tree manipulation may not be possible. Another option is to build an equality checker that makes a copy of each representation in some canonical form and then compares their representations. Decomposing the expressions into their Taylor series and comparing the series could be an option. This doesn't make equality checking as cheap as may be desirable, however.

I hope this article has offered some insight into a way of working with mathematical expressions analytically. Whether it is a good way is debatable, though even if you think this entire approach is a mistake, it's another person's mistake to learn from. While you may not need to deal with mathematical expressions regularly (or ever), some of the principles illustrated here (such as the idea of using a different internal data representation that is unseen by the end user) are more widely applicable.

TPJ