# “Oldest” bug in computer algebra system?

The goal of this question is to find an error in a computation by a computer algebra system where the 'correct' answer (complete with correct reasoning to justify the answer) can be found in the literature. Note that the system must claim to be able to perform that computation, not implementing a piece of (really old) mathematics is sad, but is a different topic.

From my knowledge of the field, there are plenty of examples of 19th century mathematics where today's computer algebra system get the wrong answer. But how far back can we go?

Let me illustrate what I mean. James Bernoulli in letters to Leibniz (circa 1697-1704) wrote that [in today's notation, where I will assume that $y$ is a function of $x$ throughout] he could not find a closed-form to $y' = y^2 + x^2$. In a letter of Nov. 15th, 1702, he wrote to Leibniz that he was however able to reduce this to a 2nd order LODE, namely $y''/y = -x^2$. Maple can find (correct) closed-forms for both of these differential equations, in terms of Bessel functions.

An example that is 'sad' but less interesting is $$r^{n+1}\int_0^{\pi}\cos(r\rho \cos (\omega))\sin(\omega)^{2n+1}d\omega$$ with $n$ assumed to be a positive integer, $r>0$ and $\rho$ real; this can be evaluated as a Bessel functions but, for example, Maple can't. Poisson published this result in a long memoir of 1823.

One could complain that (following Schloemilch, 1857) that he well knew that $$J_n(z) = \sum_{0}^{\infty} \frac{(-1)^m(z/2)^{n+2*m}}{m!(n+m)!}$$ Maple seems to think that this sum is instead $J_n(z)\frac{\Gamma(n+1)}{n!}$, which no mathematician would ever write down in this manner.

Another example which gets closer to a real bug is that Lommel in 1871 showed that the Wronskian of $J_{\nu}$ and $J_{-\nu}$ was $-2\frac{sin(\nu\pi)}{\nu z}$. Maple can compute the Wronskian, but it cannot simplify the result to $0$. This can be transformed into a bug by using the resulting expression in a context where we force the CAS to divide by it.

For a real bug, consider $$\int_{0}^{\infty} t^{-\lambda} J_{\mu}(at) J_{\nu}(bt)$$ as investigated by Weber in 1873. Maple returns an unconditional answer, which a priori looks fine. If, however, the same question is asked but with $a=b$, no answer is returned! What is going on? Well, in reality that answer is only valid for one of $0\lt a\lt b$ or $0\lt b \lt a$. But it turns out (as Watson explains lucidly on pages 398-404 of his master treatise on Bessel functions, this integral is discontinuous for $a=b$. Actually, the answer given is also problematic for $\lambda=\mu=0, \nu=1$. And for the curious, the answer given is $$\frac{2^{-\lambda}{a}^{\lambda-1-\nu}{b}^{\nu} \Gamma \left( 1/2\nu+1/2\mu-1/2\lambda+1/2 \right)} { \Gamma\left( 1/2\mu+1/2\lambda+1/2-1/2\nu\right) \Gamma \left( \nu+1 \right)} {F(1/2-1/2\mu-1/2\lambda+1/2\nu,1/2\nu+1/2\mu-1/2\lambda+1/2;\nu+1;{\frac {{b}^{2}}{{a}^{2}}})}$$

EDIT: I first asked this question when the MO community was much smaller. Now that it has grown a lot, I think it needs a second go-around. A lot of mathematicians use CASes routinely in their work, so wouldn't they be interested to know the 'age' gap between human mathematics and (trustable) CAS mathematics?

-
This is pretty much a duplicate of this thread: mathoverflow.net/questions/11517/computer-algebra-errors – Ryan Budney Feb 25 '10 at 4:34
Can you give an example to get things started? – Steve D Feb 25 '10 at 4:35
@Ryan: I don't agree - that other thread threw up excessively wide-ranging answers and so was unfocussed. This is a clearly better question. It is possible that it will throw up a clearly correct historical source, but I rather think it won't, though, and would better be community wiki. – Charles Stewart Feb 25 '10 at 11:48
Hi Charles. I suspect this question is maybe too focused. Presumably the author made a typo, for example, as I'm not awareof any 19th century computer algebra packages, as this is before the electronic computer. Perhaps this question has a simple answer. The first computer algebra package is reportedly Schoonschip (1963). Presumably it wasn't error free? – Ryan Budney Feb 25 '10 at 13:06
The only way I can see one could be surpised by that is that one very. very grossly underestimate the difficulty of simulating, say, Bernoulli and that one has never wrote any code (for anyone who's programmed anything knows bugs are essentially inevitable!) – Mariano Suárez-Alvarez Feb 25 '10 at 18:26

In Minkowski space-time one expects $$\epsilon_{ijkl}\epsilon_{i^\prime j^\prime k^\prime l^\prime}g^{ii^\prime}g^{jj^\prime}g^{kk^\prime}g^{ll^\prime}=-24, \tag{1}$$ where $\epsilon_{ijkl}$ is the Levi-Civita tensor and $g^{ij}$ represents the metric. However, if you (uncritically) calculate the l.h.s of (1) in the symbolic manipulation system FORM using FixIndex statement to assign specific values to selected diagonal elements of the Kronecker delta, which by default represents the metric, you get +24, not -24.

In fact, this is not a bug but a subtlety of inner workings of FORM and the users were warned that one can try to change the behaviour of the Kronecker delta a bit but "this is dangerous and needs, in addition to a good understanding of what is happening, good testing to make sure that what the user wants is indeed what does happen": http://www.nikhef.nl/~form/maindir/documentation/reference/online/

Interestingly, this subtlety of manipulations with the Levi-Civita tensor in FORM led to a long lasted sign error in the calculations of the pion-pole dominant term in the hadronic light-by-light scattering contribution to the muon anomalous magnetic moment, and the source of this error was discovered only when some later time calculations of the same quantity, based on the REDUCE Computer Algebra System, gave an oposite sign: http://arxiv.org/abs/hep-ph/0112102 (Comment on the sign of the pseudoscalar pole contribution to the muon g-2, by Masashi Hayakawa and Toichiro Kinoshita).

I think the following paragraph in the above cited FORM manual gives a very good advice how modern computer algebra systems should be used:

"As in the Zen saying:

To the beginning student mountains are mountains and water is water. To the advanced student mountains stop being mountains and water stops being water. To the master mountains are mountains again and water is water again.

Of course the modern master also checks that what he expects the system to do, is indeed what the system does."

-
From looking at the linked arxiv paper it seems that the problem is in that the Levi-Civita symbol is defined relative to the Kronecker delta, so that changing one changes the other? (I tried to understand your first two paragraphs and am not sure if I got it right.) If that's the case this is less of a bug in the software but a PEBCAK. – Willie Wong Oct 15 '14 at 11:21
FORM does not have an internal metric tensor. Instead it has a Kronecker delta d(mu,nu) which will produce a dot product of two vectors in p(mu)*d_(mu,nu)*q(nu) when mu and nu are summable indices. All dependence on the metric should be provided and controled by user and this is a potential source of errors. For example, if a Kronecker delta's survive in the output, they are usually interpreted as a metric tensor and in many cases this works out well. – Zurab Silagadze Oct 16 '14 at 6:22
However, the contraction of two Levi-Civita tensors will give products of Kronecker delta's and they are indeed Kronecker delta's and not metric tensors. So you will get a wrong answer if you interprete them as metric tensors. There is no real bug in the software. The problem is that to get correct results in seemingly simple manipulations with the Levi-Civita tensor the user must understand very well how FORM actually makes its calculations. – Zurab Silagadze Oct 16 '14 at 6:44

According to Wolfram Alpha and the tables in [2], $\pi(10^{10}) = 455, 052, 511$. Nevertheless, in Don Zagier's paper listed below we find that $\pi(10^{10}) = 455, 052, 512$.

Wonder whether someone has already noted this discrepancy between the sources elsewhere. Naturally, the discrepancy implies the existence of a bug in either the routines of Zagier or in WA's implementation of the prime counting function. I don't think that it's only a typo in Zagier' note because, if my memory serves me right, there are some other texts in the literature that endorse the computations of Zagier (for instance, see [1, page 7].).

References

[1] A. E. Ingham. The distribution of prime numbers. Cambridge Mathematical Library, 1934 (Reissued in 1990).

[2] H. Riesel. Prime Numbers and Computer Methods for Factorization. Birkhäuser, Second Edition, 1994.

[3] D. Zagier. The first 50 million primes. Math. Intelligencer, 0 (1977).

-
The discrepancy has been noted elsewhere, e.g., oeis.org/A006880 which says "Lehmer gave the incorrect value 455052512 for the 10th term." Unfortunately, it gives no Lehmer citation, and doesn't even say which Lehmer! It has been noted elsewhere at MO that D N Lehmer, father of D H Lehmer, was one of the last mathematicians to count 1 as a prime, so maybe there's no bug here, just two different definitions. – Gerry Myerson Sep 16 '10 at 3:13
My speculation in my previous comment was off the mark. J C Lagarias, V S Miller, and A M Odlyzko, Computing $\pi(x)$: the Meissel-Lehmer method, Math Comp 44 (1985) 537-560, say "[D H] Lehmer used an IBM 701 to calculate $\pi(10^{10})=455,052,512$ (a value later shown [by J Bohman, On the number of primes less than a given limit, BIT 12 (1972) 576-577] to be too large by 1)." This paper also notes that Meissel's calculation of $\pi(10^9)$ was too small by 56, and Bohman's value of $\pi(10^{13})$ was too small by 941! – Gerry Myerson Sep 16 '10 at 3:34
Don't think it's just a matter of definitions... All of the other entries in the table would be wrong in that case! – J. H. S. Sep 16 '10 at 5:06
Is it really that hard to go all the way up to $10^{10}$ with brute force? – Andrej Bauer Oct 15 '14 at 8:23
@AndrejBauer: Running Eratosthenes sieve up to $10^{10}$ on a modern laptop won't be a problem at all, but back in 1970s equipment was much different... – Michael Oct 15 '14 at 16:27

In mathematica, if you look at the dirichlet characters modulo 4, you don't actually get the characters.

-
What year would you associate to this? – Jacques Carette Sep 15 '10 at 16:07
I think Dirichlet introduced characters in the 1830s, in his proof of the existence of primes in arithmetic progressions. – Gerry Myerson Sep 16 '10 at 3:49

If I recall correctly from ~30 years ago, on the Apple ][ the calculation 7^2 would return 49.0001. More than 100 years ago (or even 100 years before that), mathematicians already knew that the square of an integer is an integer.

-
Hmmm. I certainly don't remember anything of the kind. Do you mean in Applesoft BASIC? If you go to <a href="calormen.com/applesoft/">this</a>; Applesoft interpreter and run the program "10 PRINT 7^2", then you will get 49. – JBorger Sep 16 '10 at 9:17
Let's try this again: calormen.com/applesoft – JBorger Sep 16 '10 at 9:18
I was good with computers when they had 6502 processors, but then it all changed... – JBorger Sep 16 '10 at 9:19
and, until you prove me wrong, -1 for sullying The Woz's name! – JBorger Sep 16 '10 at 9:23
This problem is a bit wider than that: most systems use floating point for calculations instead of integers. Even today, when most systems would use integer arithmetic for integers and floating point for non-integer reals, you can force the same "error" with an expression that explicitly leads the system into the floating point realm, for example 7.0^2.0 – Michael Oct 15 '14 at 14:26

I'm not sure I completely understand the question, but many versions of Maple give the wrong answer when counting the number of partitions of $n$, for some $n$. Obviously, mathematicians have known how to do this since at least Euler. (One could argue that mathematicians have known how to count for a very long time, indeed.)

-
This is closer. An actual answer would have quoted a paper of Euler's which contains a specific enough method that, if implemented on today's computer, would correctly and efficiently compute the number of partitions of n. The sizes where Maple returns the wrong answer are large enough that few humans would have ever computed these without mechanical help (with 'computers' like Gauss excepted). – Jacques Carette Feb 25 '10 at 20:41
What is typically attributed to Euler is the [pentagonal number theorem][1] which gives a recurrence for computing the number of partitions of n. This could certainly be implemented easily on today's computers and I suspect (though haven't tried) that it would work reasonably fast on the numbers in question. That said, I think most 'state-of-the-art' algorithms use [Rademacher's formula][2] which is more efficient and was certainly not known to Euler. [1]: en.wikipedia.org/wiki/Pentagonal_number_theorem [2]: en.wikipedia.org/wiki/Partition_(number_theory) – Jason Bandlow Feb 25 '10 at 21:09
And here is a link to Euler's work: front.math.ucdavis.edu/math.HO/0510054 – Jason Bandlow Feb 25 '10 at 21:13
@Jacques Carette: As mentioned in the link to the OEIS, you can at least turn this into an example known by 1920. Simply compute p(11269) mod 5. An identity by Ramanujan says that is 0, while Maple would say 1. – aorq Mar 28 '10 at 5:23

I don't know if you mean this but have a look here (there some bugs that seem to be quite elementary):
- http://www.walkingrandomly.com/?p=801
- http://www.walkingrandomly.com/?p=578
- http://www.walkingrandomly.com/?p=88
- ...search for "bug" on this site

Hope this helps

-
Somewhat like that, but I want the link to an early paper which showed that > 100 years ago, mathematicians already knew how to get the correct answer. – Jacques Carette Feb 25 '10 at 17:33