Laurent is a manager of mathematical software development at Waterloo Maple. He can be reached at lbernardin@ maplesoft.com.
The convergence of numeric and symbolic computation systems is leading to the emergence of analytic computation systems. These systems make it possible for scientists, engineers, and programmers to enjoy the best of both worlds -- the speed of numeric computation and the flexibility and accuracy of symbolic computing. Among the computation systems that fall into this category are MATLAB from The MathWorks (http://www .mathworks.com/), as well as symbolic computation systems such as Maple (http://www.maplesoft.com/) from Waterloo Maple (the company I work for) and Mathematica (http://www.wolfram .com/) from Wolfram Research. MATLAB, for instance, has added symbolic functionalities by bundling a stripped-down version of the Maple engine within its Symbolic Math Toolbox; while Maple has integrated numerical libraries from the Numerical Algorithms Group (http://www.nag.co.uk/) into its Linear Algebra package.
A key feature of an analytic computation system is its access to vast libraries of mathematical algorithms that perform complex tasks. These tasks range from differentiating or integrating large expressions symbolically and solving differential equations numerically or in closed form, to solving equations involving arbitrary expressions or computing a Taylor series to a high order.
As is typical of analytic computation systems, Maple combines a GUI, programming language, and large library of scientific and mathematical algorithms to perform tasks such as finding Eigenvalues of large numeric matrices using arbitrary precision, or computing exact indefinite integrals. (All in all, Maple's library includes more than 4000 functions, totaling over half a million lines of code. The Maple library is written using its own user-level programming language, making it possible for users to inspect and extend the library.) In this article, I will provide examples of the methods and data structures typically used by analytic computation systems. For the purposes of illustration, I'll refer to Maple Version 6, which runs on Windows, Macintosh, and UNIX (including Linux).
Hybrid Symbolic-Numeric Techniques
To see how arbitrary precision computations work, as well as show how combining exact and approximate computations can be used to gain a practical advantage, consider the well-known Fibonacci numbers as defined by the recurrence in Example 1(a). Using Maple, you can easily write a program to compute the nth Fibonacci number
> F:=n->F(n-1)+F(n-2); F(0):=1;F(1):=1;
> F(20);
10946
Unfortunately, while this program computes up to the 20th Fibonacci number, you can't use it for computing F100 or even F1000, because my algorithm takes exponential time. To correct this, you could write a nonrecursive program:
> F:=proc(n)
(r0,r1):=(1,1);
for i from 2 to n do
(r0,r1):=(r1,r0+r1)
end do;
r1
end;
The iterative version of the program is able to compute F10000, a number with 2090 decimal digits in less than one second on my machine. The algorithm went from exponential to linear time complexity. The disadvantage is that, given the recurrence formula by which Fibonacci numbers are defined, writing the iterative program is less natural than writing down the recursive version.
The recursive program has exponential complexity because, at every step of the recursion, it recomputes every intermediate Fibonacci number, all of which have already been computed in the previous step of the recursion. To compute F3 you have to compute F1 and F2. Then, to compute F4, you have to compute F2 and F3, both of which have already been computed before.
Because scientific calculations often involve recomputing results and rewriting a recursive program into an iterative one is often difficult, Maple uses a clever (though memory intensive) method of dealing with such calculations. In fact, you can rewrite the program for computing Fibonacci numbers as follows:
> F:=proc(n) option remember;
F(n-1)+F(n-2)
end proc;
F(0):=1; F(1):=1;
> F(300);
35957932520658356096176566517218909902 367214309267232255589801
The difference between this program and the first recursive program is the addition of the keyword option remember. Another difference is that this program, while recursive and easy to derive from the recurrence formula, is as efficient as the iterative version and runs in linear time. The key is that option remember instructs Maple to remember all calculations such that, when the program is called with an argument that it has already seen, the prior result is returned without doing any computation; see Figure 1. This useful feature is implemented by associating a hash table with the program to store previous arguments along with the result. Thus, calling a program with the same arguments a second time returns in constant time, no matter how expensive the first call was.
Again, you can now compute F10000. Can you do better? How about F100000? With numbers this size, the integer additions are no longer of constant cost; they are linear in the number of digits. Also, the recursive program is likely to run out of stack space because it involves 100,000 nested function calls. You can use the iterative program to compute F100000 in about 45 seconds. However, there is a faster way. Taking advantage of symbolic computing, you first derive a closed-form formula, see Example 1(b), for the nth Fibonacci number:
> rsolve( { F(0)=1, F(1)=1, F(n)=F(n-1)+ F(n-2) }, F(n) );
Now you just plug in any number for n and you get the nth Fibonacci number -- sort of. You can plug in any number for n, but transforming the result into an integer number involves symbolically expanding terms as in Example 1(c). This approach is computationally intensive and takes much longer than the iterative program. This seems rather disappointing -- Maple succeeded in computing a closed-form solution that is completely useless.
Here is where numeric computing comes to the rescue. First, notice that for n=100000, the second term does not contribute anything to the result because it is too small, see Example 1(b). Thus, all you have to compute is the first term. While computing powers of exact, symbolic expressions is hard, computing numeric powers is as easy as an=enlna. If you evaluate the first term of the formula in Example 1(b) using double-precision arithmetic, you get the first 15 digits of F100000 correctly. But F100000 has over 20,899 digits! Using arbitrary-precision arithmetic, it is possible to evaluate the formula using 21,000 decimal digits and to get the complete answer faster than with the recursive or the iterative program.
> round( evalf(eval(Fn,n=100000), 21000) );
4202692702995154386319005 [20860 more digits] 79669707537500
How should you estimate how many digits are needed to get the correct result? You can get a good enough approximation by computing the base 10 logarithm with low precision:
> evalf( log[10](eval(Fn,n=100000)) );
20898.62355
It is necessary to add a few guard digits to avoid roundoff errors propagating into the significant digits of the result.
Exact Computation
Now consider the recurrence in Example 2(a). To study this sequence, you would probably start by finding the fixpoints. In Maple you can achieve this as follows:
> solve( a0=111-1130/a0+3000/a02, {a0} );
{a0=5}, {a0=6}, {a0=100}
This sequence has three fixpoints: 5, 6, and 100. For any starting point, the sequence will either converge to one of these fixpoints or not converge at all. Pick a random starting point in Example 2(b). Now you can study the behavior of the sequence using these initial values. You can write a quick program to numerically compute the first 20 values of the sequence; see Listing One.
At first sight, it looks as though the sequence converges to 100 for this starting point. A second glance reveals that, for the first few iterations, the sequence approaches 6; but at the 10th term, something happens and the sequence jumps to 100. Figure 2 also illustrates the jump around the 10th term. Thus, there is reason to mistrust this result. Luckily, Maple supports arbitrary precision arithmetic, so you can redo the computation using 30 digits; see Listing Two.
Now you know for sure that the computations using 10 digits were wrong and that they did not use enough precision. It again looks as though the sequence converges to 6 (Figure 3 shows a nice smooth curve approaching 6). You didn't quite compute enough terms to reach 6, so compute the next 10 terms, as in Listing Three. But, according to Figure 4, something went wrong again. Everything seems fine at first, approaching 6; then something happens and the sequence converges to 100 (at 30 digits, we would in fact reach 100 exactly at a55). At this point, it is time to stop trying to solve this problem numerically. The sequence exhibits chaotic behavior around 6 and you would not be able to reach 6 from the given starting point using any (finite) number of digits. The more digits you use, the farther you will get; but something eventually happens and you get diverted to 100. Even though a numerical analyst once told me that since you can't reach 6 using floating-point arithmetic, the correct result is in fact 100, let's turn to exact computing.
The exact third term in the sequence is shown in Example 2(c). Start by computing the first 20 terms of the sequence using exact rational arithmetic as in Listing Four. Because the resulting fractions have quite large numerators and denominators and are thus hard to read, you display only an approximation (to 10 digits) for each term.
Finally, the sequence converges to 6, although slowly (as illustrated in Figure 5). As you compute more terms, you will get closer to 6 and will not experience any jumps due to roundoff errors. This is a nice example of how exact computing can take over where numerical computing fails.
Rapid Prototyping, Automatic Differentiation, and Code Generation
Most of Maple's mathematical knowledge is contained in a library written in its built-in programming language. The language, which feels a little like Pascal, lets you program and test code in an interpreted environment. Additionally, you have access to advanced data types such as lists, sets, arbitrary precision floating-point numbers or integers, polynomials, matrices, and so on. A garbage collector takes care of memory management. The combination of these features means that Maple can be used as a rapid prototyping environment for scientific applications.
For instance, first construct a 4×4 symmetric Toeplitz matrix with symbolic entries; see Example 3(a).
> T:=LinearAlgebra[ToeplitzMatrix]([x,y,z,t], symmetric);
Now compute the determinant of this matrix and factor it; see Examples 3(b) and (c).
> d:=LinearAlgebra[Determinant](T);
The next step is to construct a program that takes values for x, y, z, and t as input and produces the value of the determinant at this point.
> f:=unapply(d,x,y,z,t):
> f(1,2,3,4);
-20
> f(8,7,6,5);
52
How expensive is a call to this program? I can estimate this cost by counting the number of arithmetic operations that are involved.
> codegen[cost](f);
11 additions + 43 multiplications
Programs are first-class objects in Maple. This means that a program can be passed as an argument to another program that manipulates it and returns a modified program. For example, you can optimize the program with respect to the number of operations used; see Listing Five.
You can see that Maple has identified common subexpressions and reduced the number of multiplications needed to compute this determinant from 43 to 28. The resulting code can be executed as easily as the original, but will execute faster.
> g(1,2,3,4);
-20
> g(8,7,6,5);
52
Now I want to do something a little bit more complicated. I already have a program that takes four input variables and computes the determinant of the corresponding Toeplitz matrix. An algorithm called "automatic differentiation" can take this program as input, then output another program that computes the derivative with respect to one of the variables. In fact, I can instruct Maple to output a vector consisting of the derivatives with respect to all of the input variables, as in Listing Six.
At a given point, the resulting program outputs the vector in Example 3(d). For example:
> h(1,2,3,4);
56, -60, 0, -4
> h(1,2,5,8);
152, -276, 160, -28
Once my prototyping is done, I can run the code within Maple or, for faster execution, generate C or Fortran code (Listing Seven is generated C source code for the gradient program).
Conclusion
The examples presented here show that symbolic computation adds true value to the more traditional numeric methods. Analytic computation systems, with support for both high-speed numerics and the ability to efficiently manipulate symbolic data, have the potential to enlarge the class of practical problems that an engineer or scientist can solve with the help of an automated system.
DDJ
Listing One
> a[0] := 11/2; a[1] := 61/11;
> for i from 2 to 20 do
a[i] := evalf( 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]) )
end do;
a[2] := 5.59016393
a[3] := 5.63343106
a[4] := 5.67464821
a[5] := 5.71332183
a[6] := 5.74899458
a[7] := 5.77961409
a[8] := 5.77331611
a[9] := 5.17967242
a[10] := -6.8391055
a[11] := 191.5387186
a[12] := 102.8102519
a[13] := 100.1612232
a[14] := 100.0095189
a[15] := 100.0005641
a[16] := 100.0000335
a[17] := 100.0000020
a[18] := 100.0000001
a[19] := 100.0000000
a[20] := 100.0000000
Listing Two
> for i from 2 to 20 do
a[i] := evalf( 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]), 30 )
end do;
a[2] := 5.5901639344262295081967213114
a[3] := 5.6334310850439882697947214080
a[4] := 5.6746486205101509630400832965
a[5] := 5.7133290523805155490321989960
a[6] := 5.7491209197026380437051448349
a[7] := 5.7818109204856155794683379136
a[8] := 5.8113142382939957232038462063
a[9] := 5.8376565489587119615631669961
a[10] := 5.8609515225161319729296295759
a[11] := 5.8813772158414186066348580484
a[12] := 5.8991539057900653804656369740
a[13] := 5.9145249506789843093265503806
a[14] := 5.9277414077768100768697258775
a[15] := 5.9390504854613682145398481297
a[16] := 5.9486874924846266324223440920
a[17] := 5.9568707319889872021880625846
a[18] := 5.9637987220073017812321934030
a[19] := 5.9696491639651381786739608575
a[20] := 5.9745793622911322066401902108
Listing Three
> for i from 21 to 30 do
a[i] := evalf( 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]), 30 )
end do;
a[21] := 5.9787313055716646903432098926
a[22] := 5.9823017419650597836243768506
a[23] := 5.9866906078267980996399642244
a[24] := 6.0136522709040277799769525928
a[25] := 6.4232153924054698100400652873
a[26] := 12.7415627706574382005210260605
a[27] := 58.9699458767796995910794892852
a[28] := 95.8304070556576619076782493096
a[29] := 99.7392043815332393272183617503
a[30] := 99.9843246386708478051625918056
Listing Four
> for i from 2 to 100 do
a[i] := 111-1130/a[i-1]+3000/(a[i-1]*a[i-2]);
printf("a[%d] := %a\n",i,evalf(a[i]))
end do:
a[2] := 5.590163934
a[3] := 5.633431085
a[4] := 5.674648621
a[5] := 5.713329052
a[6] := 5.749120920
a[7] := 5.781810920
a[8] := 5.811314238
a[9] := 5.837656549
a[10] := 5.860951523
a[11] := 5.881377216
a[12] := 5.899153906
a[13] := 5.914524951
a[14] := 5.927741408
a[15] := 5.939050485
a[16] := 5.948687492
a[17] := 5.956870732
a[18] := 5.963798721
a[19] := 5.969649144
a[20] := 5.974579029
...
a[100] := 5.999999988
Listing Five
> g := codegen[optimize](f);
g := proc(x, y, z, t)
local t9, t18, t2, t1, t3, t6, t28;
t1 := x^2; t2 := t1^2; t3 := y^2; t6 := z^2;
t9 := t^2; t18 := t3^2; t28 := t6^2;
t2 - 3*t1*t3 - 2*t1*t6 - t1*t9 + 4*x*t3*z + 4*y*z*x*t
+ t18 - 2*t3*y*t - 2*t3*t6 + t3*t9 - 2*t6*y*t + t28
end proc
> codegen[cost](g);
7 storage + 7 assignments + 28 multiplications + 11 additions
Listing Six
> h := codegen[GRADIENT](g);
h := proc(x, y, z, t)
local df, t1, t18, t2, t28, t3, t6, t9;
t1 := x^2; t2 := t1^2; t3 := y^2; t6 := z^2;
t9 := t^2; t18 := t3^2; t28 := t6^2;
df := array(1 .. 7);
df[7] := 1; df[6] := 1; df[5] := -t1 + t3;
df[4] := 2*df[7]*t6 - 2*t1 - 2*t3 - 2*y*t;
df[3] := 2*df[6]*t3 - 3*t1 + 4*x*z - 2*y*t - 2*t6 + t9;
df[2] := 1; df[1] := 2*df[2]*t1 - 3*t3 - 2*t6 - t9;
(4*t3*z + 4*z*y*t + 2*df[1]*x,
4*z*x*t - 2*t3*t - 2*t6*t + 2*df[3]*y,
4*x*t3 + 4*x*y*t + 2*df[4]*z,
4*y*x*z - 2*t3*y - 2*t6*y + 2*df[5]*t)
end proc
Listing Seven
Toeplits matrix
> codegen[C](h,ansi);
void h(double x,double y,double z,double t,double crea_par[4])
{
double df[7];
double t1;
double t18;
double t2;
double t28;
double t3;
double t6;
double t9;
{
t1 = x*x;
t2 = t1*t1;
t3 = y*y;
t6 = z*z;
t9 = t*t;
t18 = t3*t3;
t28 = t6*t6;
df[6] = 1.0;
df[5] = 1.0;
df[4] = -t1+t3;
df[3] = 2.0*df[6]*t6-2.0*t1-2.0*t3-2.0*y*t;
df[2] = 2.0*df[5]*t3-3.0*t1+4.0*x*z-2.0*y*t-2.0*t6+t9;
df[1] = 1.0;
df[0] = 2.0*df[1]*t1-3.0*t3-2.0*t6-t9;
crea_par[0] = 4.0*t3*z+4.0*z*y*t+2.0*df[0]*x;
crea_par[1] = 4.0*z*x*t-2.0*t3*t-2.0*t6*t+2.0*df[2]*y;
crea_par[2] = 4.0*x*t3+4.0*x*y*t+2.0*df[3]*z;
crea_par[3] = 4.0*y*x*z-2.0*t3*y-2.0*y*t6+2.0*df[4]*t;
return;
}
}
