# “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

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

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

According to Wolfram Alpha and the tables in [2], $\pi(10^{10}) = 455, 052, 511$. Nevertheless, in Zagier's paper 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 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.

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

[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

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