
The Misfortunes of a Trio of Mathematicians Using Computer Algebra Systems [pdf] - GFK_of_xmaspast
http://www.ams.org/notices/201410/rnoti-p1249.pdf
======
sold
If you ask Mathematica whether x^n + y^n = z^n has solutions for n>2 and
x,y,z>0, the system will simplify it to "False", displaying knowledge of
Fermat's last theorem. This is taken from documentation (last example in
[http://reference.wolfram.com/language/ref/FullSimplify.html](http://reference.wolfram.com/language/ref/FullSimplify.html)).

I asked about x^n + 1 = z^n, which is a simpler special case with y=1. The
system no longer recognized it to be False. So the theorem was programmed as
thoughtless pattern matching. I think one day computers might become authentic
"creative" tools for mathematicians (as opposed to "computational" tools), but
Mathematica's philosophy seems to be a dead end in this regard.

~~~
taeric
To be fair... I would imagine you could catch more than a few humans with the
same trick. Specifically folk that might recognize the textbook Fermat's
theorem, and miss the special case.

(Sadly, I can offer no data point, as I probably wouldn't have recognized
either case...)

~~~
judk
y=1 isn't just a special case, it is an absolutely trivial case that is easily
shown to have only a degenerate solution.

~~~
taeric
Apologies, I did not mean special as in edge or hard. I meant simply that it
was a specific case of the same thing.

Regardless, I don't think that changes my point.

------
timmclean
I think the interesting stuff starts on page 2:

 _With the help of Mathematica, one of us found some counterexamples to our
conjectures. Fortunately, another one of us was using Maple and, when checking
those supposed counterexamples, found that they were not counterexamples at
all. After revising our algorithms from scratch, we concluded that either the
computations performed with Mathematica or the computations performed with
Maple had to be wrong. Things started to become clear when the colleague using
Mathematica also found some “counterexamples” to the above-mentioned result of
Karlin and Szego for the case in (1) and, even more dramatically, his
algorithm yielded different outputs given the same inputs. Our conclusion was
that Mathematica was computing incorrectly._

~~~
angersock
And because it's all closed-source, you'll never know _how_ it was finding the
wrong answers!

I think that the use of closed-source software in any public research or
scientific undertaking should probably be verboten.

~~~
hyperbovine
That would be Sage, but (and I say this as somebody who uses Sage wherever
possible) it's light years behind Mathematica and always will be. Do _you_
want to be the guy in charge of implementing all of the incredibly complicated
functionality in Mathematica? Or even better, writing quality documentation
for all of it, which Mathematica has? There is a lot of software out there for
which the intersection of the people who possess the technical knowledge to
write it correctly, and the people who are willing to do so for free, is
empty. (In the cases where it's not, we call those "graduate students" :)

~~~
octonion
I just used Sage to find all sets of 1 quadrillion positive integers that have
the same product and sum. Feel free to try that in Mathematica.

~~~
jessaustin
Is this a finite set? If so, what's the simplest way to show this?

~~~
octonion
I came up with my own algorithm and found several speedups, but I haven't
written an expanded article yet. I've implemented it in Julia, but Julia's
built-in factoring is (unfortunately) far slower than Sage's.

A short preliminary article:

[http://angrystatistician.blogspot.com/2014/08/learn-this-
one...](http://angrystatistician.blogspot.com/2014/08/learn-this-one-weird-
trick-and-easily.html)

Here's the Sage code:

[https://github.com/octonion/puzzles/tree/master/product-
sum](https://github.com/octonion/puzzles/tree/master/product-sum)

I'd be curious to see the speed of the corresponding Mathematica or Maple
code.

~~~
kgwgk
The Mathematica code below is roughly comparable to your vanilla version.
Using a 2013 iMac (3.2 GHz Intel Core i5) I get the following timings for
Sage/Mathematica (Mathematica gets more competitive for larger problems):

10^8 22s 65s

10^9 225s 635s

10^11 14h20 16h05

However, I found that Sage seems to calculate divisors faster than Mathematica
for large numbers. The following are the times to sum the divisors of one
million numbers starting at 10^8, 10^9, 10^10, 10^11:

Sage: 18s, 21s, 24s, 27s

Mathematica: 12s, 14s, 28s, 49s

Solutions[pSum_, pProd_] := With[{term = 1 + pSum _pProd}, With[{div =
Divisors[term]}, Map[Function[x, If[Mod[x + 1, pProd] == 0, (1 + {x, term
/x})/pProd, ## &[]]], Take[div, Ceiling[Length[div]/2]]]]]

WithNonOneValues[n_, nval_] := Last[Reap[If[nval == 0, Sow /@ Solutions[n - 2,
1], Module[{val = Array[2 &, nval], idx = 1}, While[idx <= nval, With[{pSum =
Plus @@ val + n - nval - 2, pProd = Times @@ val}, If[pSum + 2_val[[1]] <
pProd*val[[1]]^2, If[(++idx) <= nval, val[[Range[idx]]] = val[[idx]] + 1],
Map[Function[x, If[Min[x] >= val[[1]], Sow[Flatten[{x, val}], ## &[]]]],
Solutions[pSum, pProd]]; idx = 1; val[[1]] = val[[1]] + 1]]]]]]]

AllSolutions[n_] := Flatten[Map[Function[x, WithNonOneValues[n, x]], Range[0,
Log2[n]]], 2]

------
Fede_V
There is an excellent blog by the creator of Sage:

[http://sagemath.blogspot.it/2014/08/you-dont-really-think-
th...](http://sagemath.blogspot.it/2014/08/you-dont-really-think-that-sage-
has.html)

Sage is absolutely amazing, but it is very far behind the cutting edge CAS
systems (MAGMA/Mathematica) - I think he is being a bit harsh in his
assessment, but he has incredibly ambitious goals - and by those goals, he has
fallen short.

~~~
danbruc
When I first came across the Sage Wikipedia page and especially the »Software
packages contained in Sage« section [1] I immediately wondered how well that
approach works. I admittedly don't know anything about the way Sage is build
but you get the impression that they glued together many existing specialized
solution to form one comprehensive solution. And while every piece is a mature
solution and many people invested years of effort to build them, I also
experienced over and over again how hard it is to integrate heterogeneous
components. And that raises the obvious question whether it is more effective
to integrate mature but heterogeneous solutions or to build a less mature but
homogeneous solution at the added costs of reinventing the wheel.

[1]
[http://en.wikipedia.org/wiki/Sage_(mathematics_software)#Sof...](http://en.wikipedia.org/wiki/Sage_\(mathematics_software\)#Software_packages_contained_in_Sage)

~~~
williamstein
When you have a relatively tiny budget, and an absolutely massive difficult
system to build (requiring input from a wide range of specialists), there is
no possible way to make significant progress without fostering a strong
environment of community and collaboration. Integrating heterogenous
components from a wide range of high-quality specialist communities (GAP,
Pari, Maxima, etc.) as I did with Sage, brought many of these communities
together, with them work toward a common goal. We had over 60 coding spring
workshops in the last 10 years
([http://wiki.sagemath.org/#Past_and_future_workshops](http://wiki.sagemath.org/#Past_and_future_workshops)).
Also, we wrote an _enormous_ amount of new code in Sage itself, likely more
lines of code than all other math components put together.

~~~
fdej
I somehow doubt the amount of code in Sage itself is close to the sum of its
mathematical components, though Sage does contain a huge amount of new code
(most importantly a lot of high level mathematical code not found anywhere
else).

According to openhub.net
[[https://www.openhub.net/p/sage/analyses/latest/languages_sum...](https://www.openhub.net/p/sage/analyses/latest/languages_summary)],
Sage has 1.7 million lines of code. Maxima alone has 700K lines, Singular has
640K lines, SciPy has 650K lines, SymPy has 480K lines, etc.

~~~
williamstein
We haven't done "lines of code" estimates for Sage in a while, so I'm not
sure. Of course, lines of code being large isn't necessarily something to be
proud of (though there's some meaning in lines of code that have come and gone
overall). I'm very dubious about your lines-of-code estimates above. For
example, you claim "Maxima alone has 700K lines", but I just downloaded Maxima
from
[http://sourceforge.net/projects/maxima/files/](http://sourceforge.net/projects/maxima/files/),
then ran sloccount on the src directory and I get 129,584 lines of code:

    
    
        /tmp/maxima-5.34.1$ sloccount src
        ... 
        Total Physical Source Lines of Code (SLOC)                = 129,584

~~~
fdej
Yeah, I don't know what's up with that. I got the number from openhub. But
sloccount tells me that sage-6.2 only has 570K lines of code. And sloccount
says that Singular has 280K lines of code.

------
danbruc
A couple of weeks ago I noticed that Wolfram Alpha when asked for a uniquely
3-colorable graph [1] will hand you a graph that is not uniquely 3-colorable
(which is easy to verify by evaluating the presented chromatic polynomial at 3
which should yield 6 for a uniquely 3-colorable graph but actually yields 12
indicating two possible 3-colorings).

If you have a look at the Wolfram MathWorld article »Uniquely k-Colorable
Graph« [2] you see the same wrong graph with the following information.

 _It is implemented in Mathematica as GraphData[
"UniquelyThreeColorableGraph"] and Uniquely3ColorableGraph in the Mathematica
package Combinatorica` , in the latter case replacing the incorrect graph
(from Harary 1994, cover and p. 139) present in earlier version._

So this bugfix did not really fix the bug. I reported it using the feedback
form and we will see if the third attempt will get it right.

[1]
[http://www.wolframalpha.com/input/?i=Uniquely3ColorableGraph](http://www.wolframalpha.com/input/?i=Uniquely3ColorableGraph)

[2] [http://mathworld.wolfram.com/Uniquelyk-
ColorableGraph.html](http://mathworld.wolfram.com/Uniquelyk-
ColorableGraph.html)

------
fdej
Relevant:
[http://reference.wolfram.com/language/tutorial/WhyYouDoNotUs...](http://reference.wolfram.com/language/tutorial/WhyYouDoNotUsuallyNeedToKnowAboutInternals.html)

------
silentvoice
Symbolic math is not the same thing as formalized math, and doesn't carry with
it the same guarantees of correctness. That's the price one pays for having
black box functions spit out answers to hard problems in a reasonable time,
they won't spit out proofs of correctness alongside their answer.

Any symbolic system should be treated with care, but they can of course be
extremely useful. My typical use case is to have it compute very challenging
symbolic expressions for me. I then treat it as a "very plausible hypothesis"
which I then prove.

~~~
monochr
That's only true for closed source systems. You can easily debug open source
ones, look at their insides and so on.

This whole article is really just telling people to use the scientific method
in computing too: if you don't know how to reproduce a result it's not
science. Fast computation which a human can't possibly do is something that's
new from the last 30 years and most scientists still don't know how to deal
with it.

It is very sad to see that statisticians have switched to R, biologists to
python, computer scientists to the gnu tool chain, yet physics and maths seem
to have been colonized by mathematica when you have a whole set of open source
tools which are superior in every way: pari pg, gnu arbitrary precision
libraries, axiom, maxima and sage is you want everything under one roof with a
unified interface.

~~~
fddr
As a former physicist that used Mathematica heavily not that long ago, I have
to disagree. It would be great if there were open source alternatives that
were "superior in every way" but that is just not the case.

For some use cases at least, Mathematica is clearly better than anything else
I have tried.

I feel like something like ipython notebooks with the right combination of
libraries might eventually get there, but that is unfortunately still years
away.

~~~
monochr
As another physicist who used mathematica what you mean is that you were too
lazy to think about what assumptions were made implicitly in the calculations
you were performing and liked mathematica because it did them automagically
for you.

This is shown no where better in the paper than when they try to calculate
Integrate[Exp[-p _t]_ (Sinh[t])ˆ3, {t, 0, Infinity}]. In a good cas, such as
maxima, you will get no result and a ton of errors. Which is what you should
get without specifying what sort of variable p is, is it a matrix, polynomial,
group of some sort? That it's implicitly assumed to be a real number that
might turn out to be complex under some circumstances isn't a feature, it's a
bug.

~~~
fddr
I am genuinely surprise at how confidently you (quite wrongly) diagnosed my
problem.

What I did mean, was that for my use cases, I found Mathematica was superior
to the alternatives. You are welcome to think this was because I was lazy or
misguided, but it was definitely not because I liked not having to specify
domains for my variables (and I do remember having to write Assuming[p>0 && x
\in Reals, ...] and the like often, it definitely does not assume everything
is a positive real).

One thing I used Mathematica a lot for was numerical integration of ODEs (that
had been derived using the CAS part of Mathematica and had some pretty nasty
coefficients). NDSolve in my experience was just better than competition. You
can definitely get nonsense out of it, but with a modicum of care it works
incredibly well.

------
incompatible
Interesting. You may wonder how a calculation could generate different results
when run twice on the same arguments (assuming no obvious flaws such as global
variable usage). I once tracked down such a bug when calculating with bignums
in an open source program. I don't remember the details, but it was something
to do with garbage collection, perhaps not protecting part of the data from
getting collected so that it randomly became corrupt.

------
jhallenworld
Does Mathematica not have some kind of "show me the work" mode? I seem to
remember that Macsyma had this..

~~~
reikonomusha
No it doesn't. The primary reason is the way a CAS solves problems is not the
way a human does.

A secondary reason with closed systems in particular is that it may reveal
methods it uses which is a competitive advantage to other systems.

~~~
gh02t
> A secondary reason with closed systems in particular is that it may reveal
> methods it uses which is a competitive advantage to other systems.

I doubt this is the case. Mathematica at least doesn't really contain any
"magic" in the actual algorithms and uses pretty well-known methods that are
for the most part documented, if not open-source. Most of their competitive
advantage is in the design of the language itself, the interface and in their
internal virtual machine. These details aren't really things you would get any
insight into with a "show the work" mode.

Your first reason is totally right though. "Show the work" mode requires a lot
of extra stuff, because the way a person reasons through solving some thing
like an integral isn't much like how a computer solves it. To make it spit out
something really comprehensible to a person, you have to have lots of
heuristics and it isn't generally very reliable. Even automated theorem
provers (which Mathematica has some very basic functionality for) spit out
rather hard to follow or extremely long proofs. Those are much more rigorous
and less reliant on heuristics than a CAS and still don't produce "pretty"
output (AFAIK, I'm not super experienced using automated provers or proof
assistants).

Showing steps has been done before, e.g. by Macsyma like OP mentioned and by
Wolfram too (Alpha does it). It's just not that great for problems of any
substantial complexity.

~~~
htns
You could feed the proof to a (simple, verified) proof checker rather than
check it by hand. Then you only need to check the statement that was proven
was the one you wanted, without having to vet the code that spits out the
proof line by line (until you actually hit a bug and the proof fails to
check). I don't think it would be inconceivable for mathematica to implement
something like that, though perhaps it's not really their core business.

------
enupten
Makes me wonder if Timothy Daly has it right with Axiom: [http://www.axiom-
developer.org/](http://www.axiom-developer.org/)

Pity it never compiles, and has been forked into two other projects.

------
Yardlink
It struck me as odd that they avoided floating point numbers because of
impreciseness, then used some other type of number without evaluating how
precise it was. They're certainly not using typical integers with values like
10^9000. I wouldn't be surprised if those numbers are internally divided by
10^n and stored in doubles anyway - bringing them back to the same problem
they're trying to avoid.

"the determinant of bigMatrix is, approximately, 1.95124219131987·10^9762."

Approximately? How many significant digits in the exact answer? If Mathematica
promises arbitrary precision, even to 10,000 digits, they have a case. If not,
it shouldn't be any surprise.

As for getting different results on each run, I've seen this with iterative
algorithms that use random starting values. For numerically unstable problems,
that can lead to different results each time because it converges to rounding
error.

That doesn't mean Mathematica is wrong any more than getting wrong results
from doubles means its wrong. It just means it's poorly designed so the user
doesn't know about this limitation.

~~~
dalke
It's clear you aren't a mathematician. They work with "typical integers" like
10^9000 all the time, and certainly not as a trivial scaling of doubles.

There are 9763 significant digits in that calculation. They wrote it in
shorter form because the actual digits were irrelevant. The numbers they got
from Mathematica had the wrong sign and were off six orders of magnitude.
There's no reason to list more digits to show that's the case.

Mathematica promises arbitrary precision. Here's the documentation.
[http://reference.wolfram.com/language/ref/Det.html](http://reference.wolfram.com/language/ref/Det.html)
. It says "Use exact arithmetic to compute the determinant", for the
construction given in the paper. The also submitted a bug report, and got the
statement "It does appear there is a serious mistake on the determinant
operation you mentioned."

Yes, of course random algorithms can end up with different answers. The
determinate definition is not random, the documentation for Det doesn't say it
uses a random implementation, non-random methods to compute it are well known.
Again, the vendor says it's wrong - why are you disagreeing with practicing
mathematicians and the vendor? Why do you think you know enough about the
topic to be able to judge what "odd" means, in the context of this sort of
field?

------
conistonwater
Quote from the paper:

> To avoid the usual problems with floating point arithmetic (rounding,
> truncating, instability), we construct all our examples with integers.

I think this is complete overkill. I know floating-point arithmetic can be
tricky, but it's not at all random. If you have a stable algorithm and a
multiple precision floating-point library, you can use floating point numbers.
It won't necessarily give you a proof that your result (unless you spend time
to derive rigorous error bounds), but it will certainly be good enough.

I think people sometimes do slightly irrational things for fear of floating-
point arithmetic.

~~~
wglb
As someone who has done lots of work with floating point, I disagree. My fear
(or more properly caution) with floating point includes tricks like if you are
going to add up an array of floating point numbers, you need to sort them
first. The experienced floating-point programmer will know what order you need
to do this in.

The inability to represent 0.1 properly in binary floating point is not a
problem solved by increasing the precision your floating point library.

Also, there can be be surprising cases where single-precision might not be
enough. For example, if you are going to represent exchange prices for, say,
futures, you need to take into account that some of the price increments are
in pennies, some in fractional dollars. You need the calculations to be
reversible, and to have enough precision and accuracy to represent the whole
range of price increments and expected volume.

One example of a lack of appropriate paranoia about floating point numbers was
the i'm-not-very-good-at-writing-parsers' author of ASP and the result was
[http://www.exploringbinary.com/php-hangs-on-numeric-
value-2-...](http://www.exploringbinary.com/php-hangs-on-numeric-
value-2-2250738585072011e-308/) denial of service to all unpatched PHP sites.
Also bit java.

Note that interesting prime number work depends on integers and won't work
with floating point. Including factoring large numbers in cryptography.

~~~
GregBuchholz
>tricks like if you are going to add up an array of floating point numbers,
you need to sort them first.

...I think you'll want to use a priority queue, with the smallest numbers on
top, adding the first two, and then pushing the result back onto the priority
queue. Repeat until you only have one item left.

~~~
the_one_smiley
The number of people who add up arrays of floating point numbers and don't
know about Kahan summation is really incredible.

