# Computer algebra errors

In the course of doing mathematics, I make extensive use of computer-based calculations. There's one CAS that I use mostly, even though I occasionally come across out-and-out wrong answers.

After googling around a bit, I am unable to find a list of such bugs. Having such a list would help us remain skeptical and help our students become skeptical. So here's the question:

What are some mathematical bugs in computer algebra systems?

Please include a specific version of the software that has the bug. Please note that I'm not asking for bad design decisions, and I'm not asking for a discussion of the relative merits of different CAS's.

• Judging by the answers below, maybe a better question would be not "What are some bugs?" but "Which websites have the most useful/comprehensive lists of bugs?". Jan 12, 2010 at 15:36
• This is possibly some sort of record: Richard Parker told me that he once typed "isprime(2)" as his first ever query to a certain computer algebra system, and got the reply "2 is not prime". He also claimed, probably correctly, that he could find a bug in any computer algebra system within 5 minutes. Aug 9, 2010 at 14:27
• Richard's story is really surprising, because all systems I know will look up small primes in a table, so someone left off 2 from that table. It's possible, I guess, but really a silly goof. Apr 27, 2011 at 16:51
• Kevin thanks. pari has unconditional thue solver too (the second thue() in the example) and it agrees with your solutions. pari's GRH conditional solver is faster than the unconditional (in some cases the unconditional might be undoable). Do you happen to know other GRH conditional thue solver implementation?
– joro
Jul 13, 2012 at 5:19
• RE: GRH thue solver. Pari developers are investigating the problem. The thread is here: pari.math.u-bordeaux.fr/archives/pari-dev-1207/msg00008.html. Developer wrote "I do not understand where the problem (missing solution) comes from yet. It looks like a mathematical bug so far...". The pari thue code is relatively small.
– joro
Jul 16, 2012 at 5:46

I don't know any interesting bugs in symbolic algebra packages but I know a true, enlightening and entertaining story about something that looked like a bug but wasn't.$$\def\sinc{\operatorname{sinc}}$$

Define $$\sinc x = (\sin x)/x$$.

Someone found the following result in an algebra package: $$\int_0^\infty dx \sinc x = \pi/2$$

They then found the following results:

$$\int_0^\infty dx \sinc x \; \sinc (x/3)= \pi/2$$

$$\int_0^\infty dx \sinc x \; \sinc (x/3) \; \sinc (x/5)= \pi/2$$

and so on up to

$$\int_0^\infty dx \sinc x \; \sinc (x/3) \; \sinc (x/5) \; \cdots \; \sinc (x/13)= \pi/2$$

So of course when they got:

$$\int_0^\infty dx \sinc x \; \sinc (x/3) \sinc (x/5) \; \cdots \; \sinc (x/15)= \frac{467807924713440738696537864469}{935615849440640907310521750000}\pi$$

they knew they had to report the bug. The poor vendor struggled for a long time trying to fix it but eventually came to the stunning realisation that this result is correct.

These are now known as Borwein Integrals.

A video on this topic, titled "Researchers thought this was a bug," is on the 3Blue1Brown YouTube channel here.

• Which means that nobody knows Fourier analysis nowdays. Very sad and discouraging story... Jan 29, 2010 at 18:47
• @fedja I disagree. It merely shows the reasonable state of affairs that someone expert in one field, having to borrow theorems from another field, may get a surprise. Feb 16, 2010 at 19:29
• The actual person at that "poor vendor" was me. I must have spent 3 days on this problem before I figured out that Jon had tricked me. And, indeed, I am an expert in computer algebra, but do not know much Fourier analysis. But Jon's proof for why this is 'correct' is quite geometrical. Feb 17, 2010 at 4:03
• @Jacques It's wonderful to hear from the 'poor vendor'. Thanks for replying. Feb 17, 2010 at 4:19
• @Voyska No, it was not reported as a bug, just as an 'oddity' (or something like that). Jon was not mean, but playful in a devious way. He will be missed. Aug 18, 2016 at 18:18

In 1999, when I first bought an HP49G, whose major selling point was a CAS, I thought I'd try summing the harmonic series $$\sum_{n=1}^\infty \frac{1}{n}$$. I was a bit surprised to see the result 1151.8697216.

Now, I wouldn't have been too surprised if it had handled an infinite sum by just adding up "a lot" of terms until the partial sums seemed to be converging, but knowing how slowly the harmonic series grows, it wasn't plausible that it could have actually summed enough terms to get to 1151.8697216.

It turned out that it knew how to numerically compute the discrete antiderivative $$\Psi(m) := \sum_{n=1}^m \frac{1}{n} \approx \ln m + \gamma$$, and in the particular mode that it happened to be in, it would replace $$\infty$$ with the largest floating-point number it could represent, which was just under $$10^{500}$$. Indeed, $$\Psi(10^{500}) \approx 500\ln 10 + \gamma \approx 1151.8697216$$.

The story has a happy ending: after changing some flags, it returned $$+\infty$$.

• For the record, the HP49G CAS source code is actually available from the author! See www-fourier.ujf-grenoble.fr/~parisse/english.html for the download link and whatnot. Unfortunately (?), it is written in RPL, which I find a little ugly. Jan 10, 2011 at 0:05

Because the most popular systems are all commercial, they tend to guard their bug database rather closely -- making them public would seriously cut their sales. For example, for the open source project Sage (which is quite young), you can get a list of all the known bugs from this page. 1582 known issues on Feb.16th 2010 (which includes feature requests, problems with documentation, etc).

That is an order of magnitude less than the commercial systems. And it's not because it is better, it is because it is younger and smaller. It might be better, but until SAGE does a lot of analysis (about 40% of CAS bugs are there) and a fancy user interface (another 40%), it is too hard to compare.

I once ran a graduate course whose core topic was studying the fundamental disconnect between the algebraic nature of CAS and the analytic nature of the what it is mostly used for. There are issues of logic -- CASes work more or less in an intensional logic, while most of analysis is stated in a purely extensional fashion. There is no well-defined 'denotational semantics' for expressions-as-functions, which strongly contributes to the deeper bugs in CASes.

• I'm not claiming Sage is bug-free, but a lot of those bugs aren't mathematical errors -- there are plenty of compilation issues, documentation problems, etc., not to mention nearly 700 tickets classified as "enhancement" rather than "defect" -- so claiming 1582 known bugs is a little misleading. Feb 17, 2010 at 14:38
• Agreed - I have edited my response to be more precise. Feb 17, 2010 at 16:29
• That is why SAGE and other OPEN SOURCE initiatives works great here: because implementation detail are publicly available ( as well as binaries, which allows testing for free). No secret methods, no unproven improvements without peer review... Feb 25, 2010 at 15:20
• But 'peer review' has essentially established (over 15 years ago) that the fundamental design of using untyped expressions when doing symbolic calculus (i.e. computing closed-forms) is unsound. SAGE does not fix that - so the fact that it is open source and publicly available changes that how? SAGE is, by explicit design, just as bad as the closed-source systems at analysis. While I am a definite fan of open source, I do not see how, in this case, that is actually relevant. Feb 26, 2010 at 18:07
• Talking about sage, note that trac.sagemath.org/report/79 has a list of known issues which silently result in wrong results, as opposed to error messages, crashes and the likes. There are just 9 at the time of this writing. It seems to be incomplete, though.
– MvG
Jul 31, 2013 at 16:29

A quite serious error in Mathematica 7 in my opinion is that it thinks $\sqrt{x^2} =x$, not $|x|$, leading for example to 2 solutions to the following differential equation: $$y'(x) = 2 y(x) (x \sqrt{y(x)} - 1) \quad y(0) =1$$ Mathematica happily gives the following solutions: $$y(x) \rightarrow \frac{1}{(1-2 e^x +x)^2}, \quad y(x) \rightarrow \frac{1}{(1+x)^2}$$ Of course, it is a theorem that there is a unique solution to a differential equation of this type, but that doesn't mean my students hand in the wrong answer in droves...

Mathematica code: FullSimplify[DSolve[{y'[x] == 2 y[x] (x Sqrt[y[x]] - 1), y[0] == 1}, y[x], x]]

• Maple gets that one right. If the issues really is sqrt(x^2)=x, then Maple gets it right because I am the one who removed that assumption for the library in fall 1994 (Mike Monagan removed the buggy transformation from the kernel earlier that summer). Feb 26, 2010 at 16:36
• Mathematica also gets it right: Simplify[Sqrt[z^2], Element[z, Reals]] returns Abs[z], while Simplify[Sqrt[z^2]] returns Sqrt[z^2] unevaluated. The bug must be in the code of Dsolve, that somehow does not use this fact. Feb 27, 2010 at 9:10
• See groups.google.com/group/comp.soft-sys.math.mathematica/msg/… for a workaround. This is indeed a bug. Mar 1, 2010 at 14:00
• To my surprise, this is still present in Mathematica 12 Mar 21, 2020 at 13:08
• Note that $\sqrt{x^2} = |x|$ is wrong for complex $x$. Unless you tell it otherwise, Mathematica assumes $x$ is a complex variable. Nov 5, 2022 at 10:34

$2^{4^{4^4}} < 4^{4^{4^4}}$

WA: False

Update: it seems to be fixed now

• I think they have fixed it know. I tried and it says "True". Dec 17, 2013 at 2:23
• Relatedly: "does 10^(10^100) equal 10^(10^100)" returns False as of 2021. Jun 14, 2021 at 7:22
• It still returns false, as of 2022... exponents smaller than 100 return false as well
– axr
Jul 29, 2022 at 12:19
• @RavenclawPrefect WolframAlpha does not seem to be able to parse that query anymore, but trying "10^(10^100) < 10^(10^100)" and "10^(10^100) > 10^(10^100)", both return false. Jun 13, 2023 at 7:50
• @TheAmplitwist: On the other hand, the difference 10^(10^100) - 10^(10^100) apparently "cannot be determined using current methods". Jun 13, 2023 at 20:29

A friend of mine told me about his experience with Maple (version 5 or 6, I think) when dealing with matrices over $\mathbb{Q}(\sqrt{2},\sqrt{3})$. When he computed the rank and the determinant for one particular $3\times3$-matrix, he was told that the rank was 3, and the determinant was equal to zero. The answer to this paradox is, that by default, for determinants the symbolic computation methods were used for radicals, and for ranks, the floating point representations of matrix elements!

This can be thought of as either a bug or his naiveness (for not checking out how to represent elements of number fields so that floating point representations never appear), but in any case is a serious argument for treating the computer algebra software with care...

• Wait a minute, how the hell can a software written by someone with at least half a brain even try to compute rank using float computations? Feb 26, 2010 at 15:51
• I guess your assumption about at least half a brain is wrong for the team that was responsible for that routine! :-/ Feb 27, 2010 at 8:23
• I wish I knew of this example when I was teaching linear algebra! I like to emphasize that someone must write all that magic software that makes the world go round, and it may well be them. Students nod sagely. Some of them later become programmers... May 25, 2010 at 5:41
• Just this week I graded a paper where a student started out with a 3x3 matrix with 2 repeated rows, very resourcefully did some row operations including division by 3, used rounding from a calculator, and managed to find an inverse for the original matrix... I had to go back to find out what the mistake was! Jun 10, 2010 at 9:33

Sometimes a CAS cannot get the right branch of inverse trig functions when calculating integrals symbolically. See for instance: https://web.archive.org/web/20160817112947/https://pantherfile.uwm.edu/sorbello/www/classes/mathematica_badintegral.pdf

Apparently this is an unsolved problem in computer algebra.

• Correct: this is an unsolved problem. There is no theory of 'analytic' integration in closed-form, the only theory that exists is purely algebraic. And these algebraic theories all ignore branch cuts. Mar 3, 2010 at 21:36
• It does not seem to have been mentioned in the paper linked above that one reason behind the mishandling of these trigonometric integrals is in mishandling the Weierstrass substitution. I have seen the following paper: doi.acm.org/10.1145/174603.174409 but I have no idea what the state of the art is nowadays for such things. Aug 2, 2010 at 2:07
• @J.M.: as far as I know, David's paper that you refer to above is the 'most recent' on this topic. Apr 27, 2011 at 12:49
• @Jacques (sorry for the late reply): That's a bit sobering... thanks. May 8, 2011 at 21:26

Wolfram alpha is saying that the series of $$\sum_k\sin(2 k \arctan(k^2))$$ does not converge:

https://www.wolframalpha.com/input/?i=sum+sin%282+k+atan%28k%5E2%29%29

instead it converges! Seems that mathematica is only dealing with limits of functions not with limit of sequences.

Another simpler example is $$\sum_k \sin(2k \pi + 1/k^2)$$:

https://www.wolframalpha.com/input/?i=sum+sin%282k+pi+%2B+1%2Fk%5E2%29

E.

• The first example is very nice, even without considering the original question; a fun exercise in calculus. The second example is really shocking! May 16, 2011 at 0:36
• Mathematica 8.0.1 commits the same error: Sum[Sin[Pi k+1/k^2],{k,1,Infinity}] May 16, 2011 at 12:52
• WolframAlpha still claims that $\sum_{k=1}^\infty \mathrm{sin}(1/k^2)$ diverges. Jun 13, 2023 at 14:54
• But when I do this query, I get the convergence Jul 14, 2023 at 11:16

This story heard from Enrico Bombieri. I do not know if it qualifies, since it is not a CAS bug, and in addition it is second-hand. However it might be quite effective in casting doubt in the mind of your students, if that's your purpose :)

E.B. was doing some Riemann zeta zero crunching on his PC some years ago, the software he wrote seemed ok, and the next step was to run it on a mainframe to get some serious data. He was given the privilege to try it on the first Cray supercomputer. Most of the time results were nice, but every now and then he got really weird results. He and his coworkers spent some awful weeks trying to catch the bug. In the end, they cornered the problem: when the Cray divided 1 by 12 the result was a negative number...

EDIT: I double checked, it was not a Cray supercomputer but a computer based on an early iteration of the Pentium chip (I guess an IBM one), and the Pentium bug was also encountered by others of course. Sorry for the inaccuracy.

• Clearly, it calculated $1/12$ by summing the series $-1 - 2 - 3 - \cdots$. :) Mar 8, 2012 at 2:32

Here is an example in Wolfram Alpha. A student had been given the assignment of finding the limit as $n$ tends to infinity of $\frac{1}{1+\frac{(-1)^n}{log(n)}}$. He had correctly arrived at the answer 1. Now he used WA to check if he was correct. WA returned 0 (the command lim n-> inf 1/(1-(-1)^n/log(n)) ). On examining the steps, it turned out that WA had manipulated a bit and used L'Hopital on the expression $\frac{log(n)}{(-1)^n+log(n)}$.

Note that if one instead asks for the limit of $\frac{1}{1-\frac{(-1)^n}{log(n)}}$ WA correctly returns 1, using the same method one usually would.

• They must be reading this; it seems to work now. Apr 27, 2011 at 13:07

From the sage-support mailing list.

Sage 5.10 claims $$\forall a,b \in \mathbb{R}, \; \sqrt{(a+b)^2}=\sqrt{a^2}+\sqrt{b^2}$$

though it contradicts it numerically for $a=1,b= -1$.

Session:

sage: var('a,b');assume(a,'real');assume(b,'real');ex=sqrt( (a+b)^2 ) - (sqrt(a^2)+sqrt(b^2));ex
(a, b)
sqrt((a + b)^2) - sqrt(a^2) - sqrt(b^2)
sage: ex.full_simplify()
0
sage: ex.subs(a=1,b=-1)
-2


My advice is never to trust a single CAS. I only wrote one computer aided paper: I did the programming on Mathematica / Linux, and my collaborator did it on Magma / Solaris. We also made a point of not communicating while writing the programs.

• Are you more trusting of Human Algebra Systems? Jan 13, 2010 at 4:33
• Human errors when doing math tend to be very different then human errors when writing computer programs. The most famous example to how hard it is to write a correct program is binary search (see e.g. Bentley's "programming pearls", or Knuth TAOCP). It's a completely trivial algorithm, but try to code it - I promise you you'll get it wrong. Jan 13, 2010 at 7:34
• Human Algebra Systems often implements some correction algorithms based on redundancy which works with analysing analogies and symmetries. As far as I know no CAS implement something which may find error because result is not so symmetric that expected... Feb 25, 2010 at 15:17
• I personally agree wholeheartedly here. Science grows on independent confirmation, as almost any part of the supply line can fail (software, compilers, hardware). As Mike Rubinstein put it (I think): Why do we trust software to compute zeros of $L$-functions? Well, we keep tweaking the program (that is, fixing bugs) until it gives $14.1347...$ for the Riemann $\zeta$-function, and then the same for all other cases... In other words, the answer it "correct" when it matches our expected reality. :) May 25, 2010 at 6:17

This error affects all versions of Mathematica from 6 to 8. The result of a function depends on what letter is chosen for argument when calling it. In the simplest case it can be illustrated as follows:

in:

$A[\text{x_}]\text{:=}\sum _{k=0}^{x-1} x$

$A[k]$

$A[z]$

out:

$1/2 (-1 + k) k$

$z^2$

The correct answer is evidently, the later. This behavior affects not only sums but also integrals, so one have to check so that the letter user for the argument not to coincide with the index variable used for definition. In the case of recursion this becomes very difficult. The following example shows that moving a factor not dependent on the index variable out of the sum sign changes the result:

in:

A1[0,x_]:=1
A2[0,x_]:=1

A1[n_,x_]:=Sum[A1[-1 - j + n, x]*Sum[A1[j, k], {k, 0, -1 + x}], {j, 0, -1 + n}]
A2[n_,x_]:=Sum[Sum[A2[j, k]*A2[-1 - j + n, x], {k, 0, -1 + x}], {j, 0, -1 + n}]

A1[1,x]/.x->2
A1[2,x]/.x->2
A1[3,x]/.x->2

A2[1,x]/.x->2
A2[2,x]/.x->2
A2[3,x]/.x->2

A2[1,2]
A2[2,2]
A2[3,2]


out:

2
5
13

2
5
12

2
5
13

• OK, so you should not use global variables as auxiliary variables in your definitions. Now why do you call this an error of Mathematica? May 3, 2011 at 9:23
• More specifically, use A1[n_, x_] := Module[{k, j}, Sum[A1[-1 - j + n, x]*Sum[A1[j, k], {k, 0, -1 + x}], {j, 0, -1 + n}]]; A2[n_, x_] := Module[{k, j}, Sum[Sum[A2[j, k]*A2[-1 - j + n, x], {k, 0, -1 + x}], {j, 0, -1 + n}]] May 6, 2011 at 8:09

In Mathematica 7, the command

Table[DirichletCharacter[4, 2, n], {n, 0, 8}]

should return a list of values of the character with modulus 4 and index 2, evaluated at 0, 1, 2, ..., 8. Instead, it returns the decidedly non-multiplicative:

{0,1,0,-1,0,-1,0,-1,0}

Analytical Number Theoretic Functions in Mathematica are (sometimes) unreliable.

https://code.google.com/archive/p/mpmath/ is part of Sage --- https://www.sagemath.org/ --- hence you can double check Mathematica values there.

(sorry, I'm not allowed to post hyperlinks...)

This might get fixed in the future, but at the time of this writing, Wolfram Alpha gets apparently sometimes confused by logarithms of complex numbers:

Wolfram Alpha -- $$\log(1+ \frac{1}{2}i) - \log(1 - \frac{1}{2} i)$$

For reference, should the problem get fixed: it claims that $$2i = 2i\cot^{-1}(2) \approx 0.9272$$.

Curiously, the numerical approximation is correct, but the symbolic form seems to be wrong.

• This seems to be fixed now. WolframAlpha only gives $2 i \cot^{-1}(2)$ as the answer, not $2i$. Jun 13, 2023 at 7:46

A bug, this time from MATLAB. While trying to obtain:

$$\int_0^a x^2\sqrt{-x^2+ax}\,\mathrm{d}x=5a^4\pi/128$$

through:

syms x a
assume(a>=0)
int(x^2*sqrt(-x^2+a*x),x)


MATLAB get unfocused between the imaginary algebra and gets the negative value (!): $$-5a^4\pi/128$$

Original question in here. Tested in MATLAB R2014b at 2017-05-16. Edit: Present in R2018a, it was already corrected in R2018b.

This was mentioned to me by the user Umesh Shankar on Mathematics SE.

WolframAlpha incorrectly solves the following congruence equation: $$4 + 6 \times 10^n \equiv 0 \pmod{7}$$. It says that the answer is $$n \equiv 4 \pmod{7}$$. The correct answer is clearly $$n \equiv 4 \pmod{6}$$.

• This one surprised me. Jun 14, 2023 at 18:31

I found a mistake made by WolframAlpha, which I also posted on a related thread on Mathematics Educators SE.

Ask WolframAlpha to Solve 1/floor(x) < 2/x for x > 1, and it says that the solution over the reals is $$x = 99/5 \approx 19.8000$$". Rearrange the terms and ask it to Solve x < 2 floor(x) for x > 1 and it gives a different solution "$$x = 759/10 \approx 75.9000$$". Increase the lower bound on $$x$$ from $$x > 1$$ to $$x > 76$$, say, to get WolframAlpha to commit to different specific answers.

Of course, the inequality is actually true for all $$x > 1$$.

• It gives the same answer if we do not specify any condition like $x>1$ also. We must be wary of CA softwares from now. But, how do we detect these mistakes in a large and complicated numerical interation? Jul 14, 2023 at 10:51

If you are performing numerical computations, then a more likely source of error is in roundoff or over/underflow. In these cases, I wouldn't say that the CAS is necessarily in the wrong, just that you need to know the properties of the underlying algorithm and either recast your input or reimplement it in a more numerically robust way. In such cases, decent introductions to numerical analysis should give you a feel for the types of issues you need to worry about.

Of course, on the matter of symbolics, then there are no excuses for errors.

David Bailey and Jonathan Borwein said in a talk yesterday that the most recent editions of both Maple and Mathematica give the nonsensical result $$\int_0^1\int_0^1|e^{2\pi ix}+e^{2\pi iy}|\,dx\,dy=0$$ I later verified this for Maple 17, entering int(int(abs(exp(2*Pi*I*x)+exp(2*Pi*I*y)),x=0..1),y=0..1).

• Absolute values are often hard for CAS's to deal with. This is not an easy problem (unless you cheat by applying some insight). Jul 16, 2014 at 1:23
• I just asked Mathematica 10 to "Integrate[ Abs[Exp[2Pi I x]+Exp[2Pi I y]], {x,0,1}, {y,0,1}]", and after some thought it returned the expression unevaluated. Jul 17, 2014 at 1:29
• I note for the record that Maple 17 has no difficulty with evalf(Int(Int(abs(exp((2*PiI)*b)+exp((2*PiI)*c)), b = 0 .. 1), c = 0 .. 1)), although of course it doesn't return the answer in the form $4/\pi$. Jul 17, 2014 at 1:54
• Mathematica 11.3 now evaluates this integral correctly as 4/pi.
– Ben
Sep 9, 2018 at 2:05
• And Maple now gets $4/\pi$. May 17, 2019 at 16:01

http://oeis.org/A110375

A110375   Numbers n such that Maple 9.5, Maple 10, Maple 11 and Maple 12 give the wrong answers for the number of partitions of n.

There are too many to be listed on the margins of MO.

Look at the archives of the newsgroups comp.soft-sys.math.maple, comp.soft-sys.matlab, sci.math.symbolic, comp.soft-sys.math.mathematica. There you can find hundreds of bugs reported.

There is a notorious CAS bug hunter who once maintained a bug list for Maple and shows more than 5000 disturbing observations. (Press the Go! button.) Or go to MapleSoft and search Maple Primes.

• The newsgroups may have hundreds of bugs, but they are buried in the midst of 10s of thousands of posts. Many are interface issues and most mathematical issues aren't bugs so much as behavior that was unexpected by a neophyte. Jan 12, 2010 at 13:36

Here are some results where different CAS give conflicting results:

1. $$\int_{y}^{\infty} \frac{e^{-x}}{x}{d x}$$ for $$y \in \mathbb{R}$$ and $$y>0$$. Wolfram Alpha gives $$\log{y}+\Gamma(0,y)$$ and sage 4.7.1 gives $$-{\rm Ei}\left(-y\right)$$

2. For all integers $$n$$, Coq proves $$n \mod 0 \equiv 0$$ and Isabelle (Wayback Machine) proves $$n \mod 0 \equiv n$$ (The proofs are just stated in theorems, I can give the exact theorems if needed). Interesting, both proofs doesn't seem to lead to inconsistency though AFAICT they depict the usual mod.

[Added] I am a fan of sage, but this bug annoyed me.

sage 4.7.2 incorrectly reports the girth of a 7 vertex graph:

H=Graph([(0, 1), (0, 3), (0, 4), (0, 5), (1, 2), (1, 3), (1, 4), (1, 6), (2, 5), (3, 4), (5, 6)])
H.girth()
4
H.is_triangle_free()
False


sage 4.3 and 4.6.2 return correct value.

sage session in the notebook and a plot of the graph

• The mod 0 example is interesting. One could argue that $mod 0$ is "not defined" and so one has a choice. E.g. some CA systems will raise an error if you use $mod 0$. Yet the definition I know is: $n \equiv m mod k$ iff $n + k\mathbb{Z} = m + k\mathbb{Z}$ (and $a mod k$ is the least integer $r\geq 0$ such that $a=kx+r$ for some $x$). This appears in many places, e.g. in group theory $C_k$ is the cyclic group of order $k$, but often $C_0$ then denotes the infinite cyclic group. Which makes sense if you set $C_k:=\mathbb{Z}/k\mathbb{Z}$. Isabelle agrees. Anybody know a why Coq does what it does? Oct 25, 2011 at 11:35
• @Max I think Coq defines $n \mod a$ as the remainder of the Euclidean division. Applied to $0$ this gives remainder of 0. Coq thinks $\frac{0}{0}=0$ too...
– joro
Oct 25, 2011 at 12:46
• @joro Division by zero in proof assistants is a feature not a bug. Jan 1 at 11:11
• @TimothyChow Thanks. I didn't claim division by zero is bug, just show incompatible uses in different software. This means code from both them can't be mixed and this might contradict human logic.
– joro
Jan 1 at 11:23

The problem in the numerical integration example is that numerical integration in Maple is done using Int, not int. The correct command should be

evalf(Int(sin(x)^44,x=0..sqrt(44)));

which should produce consistent results (and much more quickly).

Here's one I came across just now, in Maple 12. The code

with(combinat):
F := fibonacci:
phi := (1+sqrt(5))/2:
G := k -> F(k+1)/phi^k;
limit(G(n), n=infinity);


returns 0. But from the usual explicit formula for the Fibonacci numbers, which gives $F(n) \sim \phi^n/\sqrt{5}$, the output should be $\phi/\sqrt{5}$, or $(5+\sqrt{5})/10$. Replacing the built-in Fibonacci function with one that gives the explicit formula, and running the code

F := n -> 1/sqrt(5)*(((1+sqrt(5))/2)^n-((1-sqrt(5))/2)^n);
phi := (1+sqrt(5))/2:
G := k -> F(k+1)/phi^k;
limit(G(n), n=infinity);


gives the correct answer. I've encountered things like this fairly frequently when using the built-in routine for Fibonacci numbers; presumably this routine doesn't "know" the asymptotics.

• Actually, for any function unknown function P (where fibonacci is 'unknown' in this context since it returns unevaluated for symbolic arguments), <code>limit(P(n)*exp(-n),n=infinity)</code> returns 0. Polynomially descent to 0 is not enough, that will return unevaluated. I've reported the bug. The problem is that there is an implicit assumption in the implementation that unknown functions do not "grow too fast" (or go to 0 to fast or ...), which is clearly wrong. The theory for computing limits does not deal with unknown functions at all. Mar 3, 2010 at 21:28

Just found this in Mathematica 7.0 for Mac OS X x86 (64-bit) (November 11, 2008):

x=Exp[Pi Sqrt[163] ];
N[x-Round[x] ]
N[x-Floor[x] ]
N[x-Ceiling[x] ]
N[x - Round[x], 2]
N[x - Floor[x], 2]
N[x - Ceiling[x], 2]


The functions Round, Floor, and Ceiling are the obvious functions, while "N" converts the infinite-precision expression to a floating point number (the last three lines are aimed at 2-digit precision, while the first three should be 16-digit).

The first three calculations turn up as "-480." The last three give more correct values of -$7.5*10^{-13}, 1.0, -7.5*10^{-13}$.

• You do not understand. The first three should not be 16-digit. They should do the computation in 16 digit. That is something very different. Jun 10, 2010 at 8:42
• The documentation for "N" states: N[expr,n] attempts to give a result with $n$-digit precision. N[expr,n] may internally do computations to more than $n$ digits of precision.' Both N[x-Round[x],15] and N[x-Round[x],16] give the right answers, while N[x-Round[x],$MachinePrecision] gives "-480.", and$MachinePrecision itself returns "15.9546". Jun 10, 2010 at 14:19
• Indeed N[x - Round[x], MachinePrecision] returns -480.\\ But N[x - Round[x], \$MachinePrecision] returns -7.499274028018143*10^-13\\ The decision to make N[x - Round[x], MachinePrecision] mean N[x - Round[x]] is indeed a bad one. Jun 10, 2010 at 15:37 • Very peculiar. So N[x] does its work with certain precision, while N[x,k] does its work aiming for a certain precision. I've been using Mathematica a long time, and never realized this. Jun 11, 2010 at 17:23 In the paper The Misfortunes of a Trio of Mathematicians Using Computer Algebra Systems. Can We Trust in Them?, the authors report a bug in Mathematica (which is still present in the version 10) that happens when computing determinants of matrices with large integers as entries. The strangest thing of this bug is that for some matrices, the determinant function can give different values. The Mathematica notebook which reproduces the bug is available here. • That's a great article! Oct 16, 2014 at 16:03 We found some interesting bugs in Mathematica's integration software on this thread. To wit, set integral[m_,n_] = Integrate[Log[2+Cos[2Pi x]+Cos[2Pi y]] Cos[2Pi m x] Cos[2Pi n y], {x, 0, 1}, {y, 0, 1}];  Then integral[1,1] should be$1/2-2/\pi$, but Mathematica 8.0.1 returns$1/2+2/\pi$. Values for other$m$and$n$are also wrong (see the question linked above), as can be quickly verified by replacing the "Integrate" command with "NIntegrate". Curiously, if one changes the limits of integration to {x,-1/2,1/2} and {y,-1/2,1/2}, then the correct answers appear. • I don't know how Mathematica handles this integral, but I wonder if it's somehow related to the problems mentioned in Kurt's answer in dealing with branch cuts symbolically...in this case the problem is at the point$(1/2, 1/2)\$ and moving the limits around moves this to a corner instead of the middle of the interval. Mar 7, 2012 at 20:56
• This seems to have been fixed in the current version. May 18, 2017 at 16:52

The PARI/GP Thue equations solver gives wrong results when they are conditional on GRH.

Affected are at least versions 2.5.1 (latest) and 2.4.3.

? p=x^3 - 18*x^2 + 81*x + 1;a=3^3
%1 = 27
? t=thue(thueinit(p,0),a);[#t,t] \\ conditional on GRH
%2 = [3, [[0, 3], [3, 0], [19, 2]]]
? t=thue(thueinit(p,1),a);[#t,t] \\ uncoditional
%3 = [4, [[0, 3], [3, 0], [27, 3], [19, 2]]]


Found on the pari-dev mailing list http://permalink.gmane.org/gmane.comp.mathematics.pari.devel/3629.

• FYI, Mathematica has an unconditional Thue solver built in (I had something to do with that): the input Reduce[x^3 - 18 x^2 y + 81 x y^2 + y^3 == 27, Integers]'' yields the output `(x == 0 && y == 3) || (x == 3 && y == 0) || (x == 19 && y == 2) || (x == 27 && y == 3)''. Jul 12, 2012 at 14:56
• Kevin thanks. pari has unconditional thue solver too (the second thue() in the example) and it agrees with your solutions. pari's GRH conditional solver is faster than the unconditional (in some cases the unconditional might be undoable). Do you happen to know other GRH conditional thue solver implementation?
– joro
Jul 13, 2012 at 5:16
• I'm not aware of any other GRH-conditional solver. The only place I imagine GRH being relevant is in the linear-forms-in-logarithms stage, and then it would only affect the initial bound on x. If that's correct, there really isn't any reason for the implementations to be different at all. Could the GRH bound actually come out below 27 (and therefore be missing a constant or worse)? Jul 16, 2012 at 1:12
• On the mailing list Bill says this was a bug in 2.3 which was fixed in 2.5.0. I can verify that it does not affect 2.8.0. Sep 20, 2016 at 13:44