Emerging Tools for Experimental Mathematics

Jonathan M. Borwein and Robert M. Corless

1   Introduction and ``Warm-up''

If I can give an abstract proof of something, I'm reasonably happy. But if I can get a concrete, computational proof and actually produce numbers I'm much happier. I'm rather an addict of doing things on the computer, because that gives you an explicit criterion of what's going on. I have a visual way of thinking, and I'm happy if I can see a picture of what I'm working with. -John Milnor, []

Using mostly elementary examples, we discuss the use of some recent and emerging tools for experimental mathematics. The tools discussed include so-called ``inverse symbolic computation'', using lattice reduction algorithms such as ``LLL'' and ``PSLQ'', and Sloane and Plouffe's integer sequence lookup program. We concentrate on `computer-assisted discovery' of mathematical results, but a little `computer-assisted proof' creeps in as well. This paper is not a tutorial on how lattice basis reduction algorithms such as LLL or PSLQ actually work; rather, we discuss some of the ways these tools can be used to generate conjectures, and for that, a detailed understanding of the underlying algorithms is not necessary. We do hope, however, to convey some appreciation of their power.

We begin with some warm-up examples, using the Inverse Symbolic Calculator (ISC); http:// www.cecm.sfu.ca/ MRG/ INTERFACES.html. The basic idea is simple: given the first few decimal digits of some real number, we want the ISC to guess a formula for what it ``really'' is.

For example, if we input 3.14626436994198, and click on ``simple lookup'' (the default) and ``Run'', the ISC tells us that

...
3146264369941972 = (0405) 1/abs(-sr(3)+sr(2))
Your value of 314626436994198 would be here.
3146264469611207 = (0192) (5^(1/2)+4)/(exp(1/2)+1/3)
...

K = 3.146264369941972 gave the following results:
                                                      2    4
        K satisfies the following polynomial, 1 - 10 X  + X

Now consider a second ``warm-up''. If we input the number K, computed from the infinite product

As a final warm-up, consider the following two infinite products:

The ``generalized expansions'' option guesses that there is a simple generating function for the ``egyptian fraction'' of the reciprocal of (1), namely

An `egyptian fraction' is just an ordinary rational written as a sum of reciprocals of natural numbers without repeated entries in the sum. The egyptian fraction (4) corresponding to the generating function (3) is a geometric series.

So, we have seen examples where the ISC tells us something useful, tells us something incorrect, and tells us nothing.

2  A connection between the Lambert W function and Stirling's formula for n!

The Lambert W function satisfies

Figure

There is a branch point of W at x = -1/e, where W(x) = -1. See Figure 1, which can be produced in MAPLE by the command

   > plot( [ t*exp(t), t, t=-5..1 ], -1..3, -4..1 );

As our first real example of using a new tool, we look up the sequence of denominators 1, 3, 36, 270, 4320, ¼, in []. [Sometimes denominators have nontrivial common factors with numerators. Cancellation of these common factors will make any ``guessing'' procedure more difficult.] We find the sequence immediately, and the Encyclopedia gives a reference to the delightful paper [], which does not mention W or refer to any papers on W, or indeed even use it explicitly. Thus, [] would not easily be found by a normal citation search. We find out in [] that equation (6) gives coefficients needed in Stirling's formula for n!, which begins

3  Riemann Surfaces

Tools such as MATLAB and MAPLE permit easy and accurate visualization of Riemann surfaces for elementary functions [,]. Our qualitative understanding of even extremely basic mathematical building blocks can thus be affected by mathematical software tools. See [] for more discussion of visualization in general; here we concentrate on a simple technique for visualization of Riemann surfaces, namely to make 3-d plots of Âf(z) or Áf(z).

It is necessary to prove something about this technique-namely, that it really gives us a good picture of the Riemann surface and not just a 3-d plot of the imaginary part (or the real part) of the function involved. This is pursued in more detail in [], but the key point is that given w = u + iv = f(z) = f(x + iy), then we get an accurate Riemann surface by plotting, say, (x,y,v) if and only if the missing piece of information (here, u) is completely determined once x, y, and v are given. This is simple, if not quite obvious: once we have a smooth three-dimensional surface, each point of which can be associated with a unique value (i.e. ordered pair) of the map z ® w = f(z), then we have a representation of the Riemann surface of f.

This exact association is not automatic. For example, if w = ln(z) and we plot (x,y,u), then we do not get a picture of the Riemann surface for logarithm, because the branch of v = Á(w) = arg (z) is not determined from u = ln(x²+y²)/2, x, and y. If we plot (x,y,v), of course, we do recover the classical picture of the Riemann  surface for ln(z), because given x, y, and v we can easily find u.

The following short piece of MAPLE code shows how to graph the Riemann  surface for the Lambert W function. We urge you to try the following computation, because the dynamic coloured picture you get is much more easily understood than the static black-and-white image in Figure . We also urge you to try your hand at your own functions; many others are graphed in [] and [].

   > plot3d([x,y,v], u=-6..1,
            v=-5..5, axes=FRAME, 
            orientation=[-110,73], 
            labels=["x","y","v "], 
            style=PATCHNOGRID,
            colour=u,
            view=[-1..1,-1..1,-5..5],
            grid=[50,50]);

3.1  1-1 correspondence proof

If v ¹ 0, and moreover ycosv - xsinv ¹ 0, then dividing the real part by the imaginary part gives u in terms of x, y, and v:

This is precisely what we observe in the graph: two sheets intersect only if -1/e £ x < 0 (note that the colours are different and hence the corresponding sheets on the Riemann surface do not ``really'' intersect), and all sheets have a singularity at the origin, except the central one, which contains v = 0. This is as good a representation of the Riemann surface for the Lambert W function as can be produced in three dimensions.

However, Figure 2 is nowhere near as intelligible as the live MAPLE plot. On a PC, the use of OpenGL by MAPLE allows the plot to be rotated by direct mouse control. This helps to give a good sense of what the surface is really like, in three dimensions.

4  Dynamical systems, numerical analysis, and formal power series

It turns out to be useful to rescale y and t so that v = hy and t = ht, giving

The point of view taken in [] is that of backward error analysis. That is, instead of asking for the difference between v(n) and v_n, we ask instead if there is another differential equation, say

We do this not so we can improve the behaviour of Euler's method for this problem, but rather so that we may understand what Euler's method has done to the problem; for by understanding the function B(w) we will learn something about Euler's method, by comparing (9) to (11).

It turns out that we can use the method of modified equations [] to find as many terms of the Taylor series for B(w) as we desire. When we compute the modified equation for (7) to (say) fifth order, we get

We find in that paper that if

The Julia set itself approaches the origin arbitrarily closely. That is, there are poles arbitrarily close to the point of expansion of the series given for B. Thus the series (13) diverges-but, nonetheless, it can be used to evaluate B(w) for w close enough to zero, using MAPLE's built-in sequence acceleration techniques. This is precisely where the convergent infinite products are slow, and hence the series is useful. See [] for details.

But more to the point, in [], G. Labelle completely solves the problem of interpolating discrete dynamical systems with continuous dynamical systems, in the domain of formal power series. The mathematical language, however, is quite different from that used in the numerical analysis world. As an example, in [] the `modified equation' is termed an `infinitesimal generator' for the discrete dynamical system. Therefore, simple subject searches might not find []. Indeed, a combinatorics journal seems an unlikely place to find the solution of a problem in the numerical analysis of dynamical systems, but the Encyclopedia of Integer Sequences provides a way to search the `knowledge database' that is keyed on the examples, or the concrete results, of papers-not the jargon. This, if you like, is a new kind of search tool.

5  An integer-relation example

Such a conjecture might arise naturally from consideration of

To be explicit, we conjecture that

We give a short overview of using the LLL algorithm to find integer relations. Suppose that we have a finite set B of n-dimensional linearly independent vectors with rational entries. We call the set

To proceed in MAPLE, we choose a large constant C and form the following matrix, and use the lattice reduction subroutine.

We work to 30 digits for this example. In general, one has to experiment to find how many digits one needs.

We compute an approximation for the value V that we wish to identify, and approximations of the quantities that we wish to relate to V.

We now choose a large constant C. We use the size of C to penalize vectors that do not combine to zero.

We construct the rows of the matrix that we need, as follows.

   >  for i to 6 do 
   >    row.i := [ seq(0, j=1..7) ]: 
   >    row.i[i] := 1:
   >    row.i[7] := C* lastcol[i]:
   >  od:

Now we call the ``lattice'' routine to compute a short basis for the set generated by these rows.

All of these new basis vectors are of the form

Issuing the following MAPLE command lends credence to our suspicion.

There is a simpler Web-based implementation, which uses the ``EZface" to emulate a more comprehensive GNU MP implementation of this method. Go to http:// www.cecm.sfu.ca/ MRG/ INTERFACES.html, click on ``EZface", and type in the following:

lindep([1., 
9.86960440108935861883449099988, 
97.4090910340024372364403326888, 
961.389193575304437030219443653, 
9488.53101607057400712857550392, 
-16.6994737192290704961872434007])

0, -120., -140., 15., 0, -24.

Issuing the command ``lindep" calls a subroutine that looks for short integer linear dependencies among the given vector of numbers. Again its results are to be considered as ``possible relations", to be proved later.

The LLL approach is very powerful and, if used with imagination, offers rich possibilities for discovery.

One danger, indeed a `mathematical sin', is indicated by the following quotes:

``Number theory is full of things that everyone knows are true, but no one is able to prove.'' -Peter Borwein

``I have absolutely no interest in proving things that I know are true.'' -Henry Abarbanel, at a dynamical systems workshop at Penn State in 1994.

We discuss `redemption' of sins such as this in the final section.

Numerical instability in the LLL algorithm may cause difficulty, as well. Here we have simply worked to ``enough'' digits to mitigate its effects-that is, we are trying to buy more accuracy by paying for more precision. This is often expensive, and PSLQ, discussed in Section , is better.

6  How solvable is `solvable'?

RMC-removed reference to CL. Is is necessary? Should we put it back? Choosing q = exp(-2pÖ[(N/3)]) for N Î \mathbb Q, it can be shown that s_N : = c/a is an algebraic number expressible by radicals; see []. If N is a positive integer, then s_N is called the Nth cubic singular value. What can we discover computationally about s_N? For example, can we determine radical formulae for the higher order cubic singular values?

The following observations help the efficiency of the computations. It is known that

6.1  A useful transformation

The PSLQ algorithm (see []) and the LLL algorithm can both be used to find integer relations (and hence minimal polynomials for an algebraic number a, by looking for an integer relation among 1, a, ¼, a^m). However, PSLQ can also produce negative results. If PSLQ fails to find an integer relation, then one can usually say that there is no such relation with coefficients less than a computable bound, effectively proving that there is no ``simple'' relation of the guessed form.

The authors of [] used these ideas to `decode' the numerical values of s<sub>N</sub> into radical form, up to N = 100, and some values beyond, such as N = 110 and 154. They used a variety of strategies to verify the results; some ingenuity was necessary in order to extract the radicals. For N < 53, they computed P<sub>N</sub>, the minimal polynomial for G<sub>N</sub> or g<sub>N</sub>; they then tried factoring P<sub>N</sub> over different quadratic number fields until they got a factor of degree 4 or less, which they solved in radicals. This approach failed at N = 53, where they first had to use a special MAPLE program for finding a radical for any solvable quintic. (See ftp:// calfor.lip6.fr:/ pub/ softwares/ Maple/ quinticV2.gz.) The radical returned for N = 53 has over 7500 symbols in it. Kevin Hare at the CECM refined it to an equivalent but simpler radical with `only' 860 symbols. MAPLE was able to verify symbolically that this simpler radical solved P₅₃. In general, determining that a symbolic equation is indeed zero is, in certain classes of expressions, computationally undecidable [].

Indeed, the point of this whole exercise was to determine how good both MAPLE's symbolic tools and PSLQ's numerical ones were on ``grand challenge'' examples. Experience with exercises such as this have led to improvements in both tools.

If symbolic verification is not possible, reassurance that the results are correct can often be obtained by using Klein's absolute invariant []

7  Final Vignettes

Integer relation algorithms have already helped to discover many new results. We list a few of these, again taken from []. The number of such results continues to climb. We have to tell the algorithms what kind of relationship to look for, but, given that, the algorithms allow previously impossible jumps.

7.1  Zeta value series

There is a similar but more complicated formula for z(5), due to Koecher, that does suggest generalization, however. Borwein and Bradley used an LLL algorithm to determine the new coefficients []. They found that

7.2  Independent computation of digits of p

The following shows a temporary increase in complexity. This phenomenon is called ``intermediate expression swell''.

Looking at those conjugate radicals in the denominators suggests expansion-this step is natural but not automatic.

As an aside, equation (20) is an interesting identity itself. In the notation of [], it implies (once the proof is completed) that

The next step the mechanical proof of the Bailey-Borwein-Plouffe p formula simplifies the hypergeometric functions:

The next step is not necessary, but it slows down the computation so we can see that many of the terms in the above formula simply cancel.

Now, finally, our answer is plain:

A somewhat more efficient version of (19) was discovered by Fabrice Bellard. This has led Colin Percival, an undergraduate student at Simon Fraser, to design an ingenious parallel internet computation of staggeringly high order hexadecimal digits of p. Details may be found at http:// www.cecm.sfu.ca/ projects/ pihex/: the five trillionth bit of p is `0'.

7.3  Fast series for the Catalan constant

Ramanujan discovered the following series for G, which converges much more quickly than (21) []:

The general formula for Bradley's family of series is proved using certain identities among log tangent integrals. For example, (23) is proved using

David Broadhurst has, in his pursuit of new insights for theoretical physics, computationally probed more of these constants []. Based on an extraordinary blend of intuition, methodical use of PSLQ, and computer-assisted proofs, he was led to remarkable binary identities for polylogarithmic constants such as z(3), z(5), and Catalan's constant. His formula for Catalan's constant is:

Thus, digits of both G and p may be computed in the same fashion, and we might hope that the formula sheds some light on the normality of Catalan's constant. [Recall that a number is `normal' if its digits occur with equal frequency.]

8  Sin, Redemption and Cautionary Tales

Jon, eliminating this example, which I still find somewhat weak (especially when the point about proof has already been made), would also eliminate a reference and save some much-needed space. As a final demonstration, consider the power series

If L(x) and R(x) denote the left-hand and the right-hand sides of (24), respectively, then computer manipulations (under the assumption 0 < x < 1) show that dL/dx = dR/dx: mechanically differentiating both sides and using ``simplify'' reduces the difference between the two to zero. Observing that L(0) = R(0) = 0 completes a proof of (24).

8.1  Knowing `the answer' might limit us

We are all familiar with examples of the value of `doing things ourselves'. It is now trivial in most computer algebra systems (CAS) to compute very large values of the partition function with little or no thought, directly from the generating function

This really will require attention: for example, the authors of [] report in their conclusions that had they been aware of the answers in the Encyclopedia, they might not have bothered to prove what they did-and their results went beyond those in the Encyclopedia!

9  Concluding Remarks

It's not over yet. The merging of text and tools that can be anticipated over the next few years will make an enormous difference-we can expect greater insight while reading mathematical materials, and easier access to yet more powerful tools-but we make no detailed predictions, because the most significant, qualitative, changes to the work environment are by their nature unexpected.

Acknowledgements

A   The Lambert W function in brief

Here are some of the elementary properties of W.

1. On 0 £ x < ¥ there is one real-valued branch W(x) ³ 0 (see Figure 1). On -1/e < x < 0 there are two real-valued branches. We call the branch that has W(0) = 0 the principal branch. On this branch, it is easy to see that W(e) = W(1·e¹) = 1.
2. The derivative of W can be found by implicit differentiation to be

   d dx W(x) = 1 (1+W(x))eW(x) = W(x) (1+W(x))x

   where the second formula follows on using expW(x) = x/W(x), and holds if x\not = 0. We may use the first formula to find the value of the derivative at x = 0, and we see the singularity is just a removable one.
3. The function y = W( expz) satisfies

   y + logy = z .

   This function appears, for example, in convex optimization. Consider the convex conjugate, f^*(s) = sup_r rs-f(r), of the function f (r) = rln(r/(1-r))-r. Calculation shows that f^*(s) is just W(exps).

4. W(x) has a Taylor series about x = 0 with rational coefficients. Similarly, W(expz) has a Taylor series with rational coefficients about z = 1. MAPLE computes the first few terms to be

   W(ez) : = 1 + 1 2  (z - 1) + 1 16  (z - 1)2 - 1 192  (z - 1)3 - 1 3072  (z - 1)4 + 13 61440  (z - 1)5- 47 1474560  (z - 1)6 - 73 41287680  (z - 1)7 + 2447 1321205760  (z - 1)8 - 16811 47563407360 (z - 1)9 - 15551 1902536294400  (z - 1)10 + O((z - 1)11)

   Here is an exact formula for the coefficients of the nth derivative of W(expz), in terms of second-order Eulerian numbers án \atopwithdelims < > k ñ []. This formula comes from the following expression for the nth derivative of W(expz), which is stated in []. Once the answer is known, the proof is an easy induction, which we leave for the reader.

   The derivatives of W(expz) are

   dn dzn W(ez) = qn(W(ez)) ( 1 + W(ez) )2n-1  ,

   (25)

   where q_n(w) is a polynomial of degree n satisfying the recurrence relation

   qn+1(w) = -(2n-1)wqn(w) + (w+w2)qn¢(w) ,        n > 1

   (26)

   and having the explicit expression

   qn(w) = n-1 å k = 0  á n-1\atopwithdelims < > k ñ (-1)k wk+1 .

   (27)

   If n = 1 we have q₁(w) = w, and it is convenient to put q₀(w) = w/(1+w); this isn't a polynomial, but it makes things work out right. This means that our series for W(expz) about z = 1 is just

   W(ez) = å n ³ 0  qn(1) n! 22n-1 (z-1)n  .

   (28)