Algebraic Programming in Mathematica

Andrzej Kozlowski

“Algebraic Programming” means programming based on the use of algebraic concepts, such as commutativity, associativity, and distributivity. It is a method of programming for which Mathematica is particularly well suited. We give several examples of algebraic programming in Mathematica, culminating with a solution of a problem that was posted to the MathGroup email list.

I first came across the expression “the algebraic method in programming” (or algebraic programming for short) in the context of Mathematica in the book by Ilan Vardi []. By algebraic programming, Vardi meant solving programming problems by using concepts of abstract algebra, such as distributivity, commutativity, and associativity. Vardi gives some examples of algebraic programming. One of these is a solution to the following problem.

Vardi makes the observation that the answer to the problem is contained in the following algebraic formula.

Expand[s]

This formula depends only on the algebraic property of distributivity of multiplication and the existence of a unit. Indeed, Mathematica makes it possible to solve the problem directly by converting the last formula.

Variables /@ (Expand[s] /. Plus -> List) // Sort

We can write a function that will apply this method to any set of symbols. Note, however, that this program will not work if we use a list of numbers!

subSets[{a, b, c, d}]

This attractive method can be used to produce unexpected solutions to surprisingly many problems. However, this is not actually what Vardi meant by “algebraic programming.” In fact one might call this approach “naive algebraic programming.” It is naive because it makes use of specific algebraic functions rather than the general algebraic principles on which they are built. Here is a “sophisticated” algebraic program due to Vardi that achieves the same result.

subSetsV[{a, b, c, d}]

Not surprisingly, the second program not only works with symbols and numbers, but is also considerably faster than the first.

subSets[test]; // Timing

subSetsV[test]; // Timing

The main reason is that the second program is more “abstract,” in the sense that it makes fewer assumptions and has to perform fewer checks. Essentially it applies the principle of “distributivity” in the widest possible sense. In general, algebraic computations using functions like Times, Expand, and so forth involve extensive sorting of intermediate expressions before they are brought to a canonical form. From the point of view of efficiency, it is much better to operate with lists because no evaluation or canonicalization ever happens, and lists can contain anything.

. Coefficients of Polynomials and Series

Probably the best known use of algebraic programming is in solving combinatorial problems, where many solutions can be expressed in terms of so-called “generating functions.” Many examples are considered in []. Here I shall only briefly mention a problem that was sent to MathGroup in July 2001.

Here is a solution that uses algebraic programming.

Here is an example.

kpartitions[24, 5, {1, 2, 3, 4, 7, 11, 17, 19}]

The proof of the correctness of this solution is left to the reader. The above solution, while elegant, suffers from the usual defects of naive algebraic programming; it is fairly slow. Match faster solutions can be found in the MathGroup archives.

. Attributes in Mathematica

Let us look at another example. Suppose we are given a list of n objects and want to construct all permutations of k objects from them. Again, we shall look for an algebraic program to realize this task. Let us consider an example where n = 4, k = 3, and the objects are {a, b, c, d}. Then, as before, we see that we would get an answer by taking the terms of the expansion of

(a + b + c + d) (a + b + c + d) (a + b + c + 
          d)

except that this time we do not want the multiplication to be commutative. Mathematica has the following function.

?NonCommutativeMultiply

However it knows very few rules so if we are lazy we might just clear the Orderless attribute from Times and do something like this.

Perms[{a, b, c, d}, 3]

As before, this will not work for numbers instead of symbols. However, we could easily deal with this problem by mapping Hold on the list and then applying ReleaseHold at the end). Another problem is that we temporarily cleared the Orderless attribute of Times, which is generally not considered a very good idea and may produce unexpected results. In this case there were no such surprises. However, a general algebraic version of the code is both simpler and much more efficient.

Perms1[{a, b, c, d}, 3]

Here we see a characteristic difference between Mathematica and specialized programs for dealing with abstract algebraic structures like rings, groups, and so on. Such programs tend to be object oriented and contain structures like groups or rings as some form of data types. Mathematica does not have objects or types that correspond to groups or rings, but contains functions and attributes that let one implement distributivity, commutativity, associativity, and so forth in a functional way. As a bonus, we can find uses for them in more general programming, which are not strictly concerned with abstract algebraic structures.

As we have seen, distributivity is implemented in Mathematica through the function Distribute. How about the other fundamental algebraic notions, commutativity and associativity? These are essentially implemented as a so-called “Attributes” of functions Orderless and Flat. Understanding these attributes (as well as OneIdentity, which does not have a simple algebraic motivation) is the foundation of algebraic programming.

?Orderless

The attribute Orderless in a function is usually said to correspond to the mathematical notion of commutativity. This is indeed true, but one has to be careful with the interpretation of “correspond.” Many functions that are “mathematically” commutative do not have the attribute Orderless. Such is the case with the logical function Or. From the mathematical viewpoint the function is commutative.

ImpliesQ[Or[A, B], Or[B, A]]

However, Orderless is not an attribute of Or.

Attributes[Or]

We can see the effect this has on pattern matching.

a || b /. _ || a -> c

b || a /. _ || a -> c

The reason is not mathematical but has to do with performance and the way logical functions are used in programming, where one actually may wish to exploit the order of evaluation. The arguments of Or[p,q,r,s,t,...] are evaluated in order, and True will be returned the moment one of them is found to be true. This may be a substantial timesaver. Also, consider the following example, which I owe to Daniel Lichtblau.

We can see that the following two evaluations will produce different behavior.

cond1 /. a -> {x}
cond2 /. a -> {x}

In other words: we have to pay attention to order when pattern matching even though || is in fact a commutative operation. This can be exploited in programming both to increase efficiency and to avoid generating unnecessary errors.

On the other hand, it is precisely the effect of the “algebraic” attributes Orderless, Flat, and OneIdentity on pattern matching that makes algebraic programming a surprisingly powerful technique.

Let us then consider these remaining algebraic attributes, Flat and OneIdentity, which together with the function Distribute and the Orderless attribute, play the central roles in algebraic programming.

The attribute Flat, corresponds to the mathematical property of associativity. Of course, as usual there is more to it than just that.

?Flat

In order to understand the effect of the algebraic attributes on pattern matching we shall use the Mathematica function ReplaceList. Indeed, this function, together with Distribute and the algebraic attributes, is the most useful Mathematica command in algebraic programming.

?ReplaceList

Let us try using ReplaceList to understand the effect of the Flat attribute on pattern matching.

We start by taking a function f with no attributes and consider the following replacement rule.

ReplaceList[f[a, b, c], f[a_, b_] -> g[a, b]]

Not surprisingly there was no match. Now give f the attribute Flat.

Let us repeat the previous attempt.

ReplaceList[f[a, b, c], f[a_, b_] -> g[a, b]]

Note what happened. The function f[a, b, c] was now rewritten using “associativity” in two ways, as f[f[a], f[b, c]] and f[f[a, b], f[c]]. The rule was then applied but only at the top level. Now suppose we give f the attribute OneIdentity in addition to Flat.

ReplaceList[f[a, b, c], f[a_, b_] -> g[a, b]]

You see that all expressions of the type f[x] were replaced just by x. Of course this happens only when pattern matching. If you enter

f[x]

you still get just f[x], not x.

Indeed:

?OneIdentity

The useful function ReplaceAllList works like ReplaceList but also works on all subexpressions.

Let us now look at another example of algebraic programming. Of course, the most natural use of algebraic programming is solving algebra problems. So let us first consider an associativity problem.

We use a function f, which we assign the attributes Flat and OneIdentity.

FixedPoint[
      Union[Flatten[ReplaceAllList[#, f[a_, b_] -> g[a, b]]]] &, 
      f[a, b, c, d]]

Union[% /. f -> g]

This sort of thing is in fact quite useful in various areas of mathematics. Later we shall give just one application of exactly this idea to solving a rather well-known puzzle that is frequently asked on symbolic algebra email lists.

. The 4 Fours problem

The following problem was sent to the MathGroup email list in June 2002.

The problem as formulated is rather unclear because of the vagueness of the phrase “rational operator.” The question allowed not only binary operations but also unary operations, such as taking the square root or factorials. This of course means that there is an infinite number of possible solutions unless one imposes additional restrictions. However, we shall first simplify the problem and consider only binary operations: addition, subtraction, multiplication, division, and exponentiation. Let us try to solve the following simpler version of the problem.

This is very easy to implement in Mathematica using algebraic programming. In fact, we shall use exactly the same method as in the last example.

The idea is very simple. First we find all ways of writing expressions like (Subscript[A, 0])[(Subscript[A, 1])[4, 4], (Subscript[A, 2])[4, 4]], where the (Subscript[A, i]) are the binary operations Times, Plus, Divide, Subtract, and Power. Actually in Mathematica, Times and Plus are n-ary operations for any n ≥ 1, but Divide, Subtract, and Power are not. To handle them all simultaneously, we consider them to be binary operations.

It is clear that the number of different expressions that can be made up of four arguments and binary operations (Subscript[A, 0]), (Subscript[A, 1]),  ...  (which may be different or the same) can be described as follows. Let A be an associative binary operation. Then expressions such as A[x, A[y, A[z, w]]] and A[A[x, y], A[z, w]] are equal. Clearly the number of such expressions that are equal by virtue of the associativity of A (which we found in the previous section) is equal to the number of distinct ways of applying three binary functions to four arguments. In fact we can obtain all the expressions of the latter type from those of the former as follows.

First we repeat what we did above, that is, define an associative operation f.

Next, just as above, we find all the ways of “associating” four elements.

ls = 
      Union[FixedPoint[
            Union[Flatten[ReplaceAllList[#, f[a_, b_] -> g[a, b]]]] &, 
            f[a, b, c, d]] /. f -> g]

Then we remove all duplicates and expressions that involve the application of g to three elements.

ls = 
      Union[Flatten[
          DeleteCases[
            ls, _?(Not[
                    FreeQ[#, g[x__ /; Length[{x}] ≥ 3]]] &), 
            1]]]

Finally we replace g by (Subscript[A, 2]), (Subscript[A, 1]),(Subscript[A, 0]), in this order, and make all the elements {a, b, c, d} equal to 4.

ls = (i = 0; 
        MapAll[If[# === g, # /. g -> Subscript[A, Mod[++i, 3]], #] &, ls, 
            Heads -> True] /. 
          Thread[Rule[{a, b, c, d}, 4]])

(A different approach to the same problem can be found in exercise 13 of Chapter 6 of []).

We shall now simply substitute the five binary operations Plus, Subtract, Times, Divide, and Power for (Subscript[A, 0]), (Subscript[A, 1]),(Subscript[A, 2]) in all possible ways. First we need to redefine Divide and Power. Basically we do this to avoid error messages caused by division by 0. For later use, we also place a restriction on the maximum and minimum size of the exponent to avoid unnecessary computations of very large and very small numbers.

Next we create a set of rules that we shall use to create all possible ways of replacing the (Subscript[A, i]) with the binary operations.

Here is what the list we created looks like.

Take[rules, 5]

In fact, since an actual formula never involves more than three binary operations, we remove the rules involving unnecessary As.

Take[rules, 5]

Let us now find all possible solutions of the “easy” 4 fours problem.

solutions = 
      Union[Flatten[
          Union[Union[
                  Select[Flatten[Union[# /. rules]], 
                    Element[#, Integers] && 
                        0 ≤ # ≤ 99 &]] & /@ 
              ls]]]

Length[solutions]

Finally let us find all possible ways of obtaining a given number, for example, 7.

Select[patterns, 
        ReleaseHold[#] == 7 &] /. {divide -> Divide, 
        power -> Power}

This completely solves the “easy” 4 fours problem. We can solve the analogous 5 fives problem and so on in exactly the same way. In the final section we take a look at the “hard” 4 fours problem presented in the original question.

. The Hard 4 Fours Problem

The hard version of the problem is to find all possible representations of all numbers between 1 and 99 using any “reasonable” operations, including unary ones (“functions”) like Sqrt or Factorial. Of course there is now an infinite number of possible ways of trying to write out a solution, even with a finite number of operations. We shall therefore concentrate only on representations using Sqrt and Factorial in addition to the usual arithmetic operations, and only try to find some additional solutions. (If the additional operations were binary [e.g., Max or Min] then we could use exactly the same method as before to obtain a complete solution.) We first introduce modified definitions of the additional operations.

We create a new set of rules that transform our unary functions (Subscript[B, i]) into factorial, sqrt, or Identity to account for the cases where no unary operation is used.

Take[rules1, 5]

We see that there are rather many identical rules. We can reduce their number.

Take[rules1, 5]

The most obvious expressions to evaluate are obtained from the earlier ones by replacing all fours by (Subscript[B, i])[4].

Here are all the solutions we get in this way.

extendedsols = 
      Union[Flatten[Union[
                Select[Flatten[Union[Union[expr[#] /. rules] /. rules1]], 
                  Element[#, Integers] && 
                      0 ≤ # ≤ 99 &]] & /@ 
            Range[5]]]

Length[%]

An interesting case is the number 31, which was not among the solutions we got earlier. Let us see how it is represented.

Select[Flatten[extendedpatterns], 
          ReleaseHold[#] == 31 &] /. {divide -> Divide, 
          power -> Power, sqrt -> Sqrt, 
          factorial -> Factorial} /. 
      HoldPattern[Identity[x_]] -> x

Of course there are in principle infinitely many ways to insert more Bs into different slots, and it is not clear how many more numbers can be obtained in this way. Here is just one example. Consider the pattern (Subscript[A, o])[(Subscript[A, 2])[(Subscript[B, 2])[(Subscript[A, 1])[(Subscript[B, 1])[4], 4]], (Subscript[B, 3])[4]], 
        (Subscript[B, 4])[4]].

extrasols = 
      Union[Select[
          Flatten[Union[
              Union[(Subscript[A, 0])[(Subscript[A, 2])[(Subscript[B, 2])[(Subscript[A, 1])[(Subscript[B, 1])[4], 4]], (Subscript[B, 3])[4]], (Subscript[B, 4])[4]] /. 
                    rules] /. rules1]], 
          Element[#, Integers] && 
              0 ≤ # ≤ 99 &]]

We see that it produced just one new number.

Complement[extrasols, extendedsols]

However, we created several new representations for the solutions already found. Here is another representation of 31.

Select[extrapatterns, 
          ReleaseHold[#] == 31 &] /. {divide -> Divide, 
          power -> Power, sqrt -> Sqrt, 
          factorial -> Factorial} /. 
      HoldPattern[Identity[x_]] -> x

Here are all the solutions we have found so far.

solutions = Union[extendedsols, extrasols]

Here are the remaining outstanding numbers.

Complement[Range[99], solutions]

See if you can find interesting representations of these, possibly using some other operations.

Acknowledgments

I would like thank Daniel Lichtblau and Michael Trott [] for useful comments and suggestions.

References

[]  Ilan Vardi, Computational Recreations in Mathematica, Reading, MA: Addison-Wessley, 1991.

[]  R. L. Graham, D. E. Knuth, and O. Patashnik, Concrete Mathematics, Reading, MA: Addison-Wesley, 1989.

[]  Michael Trott, The Mathematica Guidebook to Programming, Berlin: Springer-Verlag, 2003.