Posted by Dave Bacon on December 29, 2007 Miguel Pais points to an interesting behavior of Mathematica, where he plots the function which is the square of the square root of x. Now, if the domain of x is taken to be complex numbers, Mathematica’s behavior seems to me to be fine. But can anyone explain this behavior
as anything other than a bug?
Update: Oops. That wasn’t the one I was trying to paste. See what happens when I disconnect from the intertubes for a few days. How about this one:
In Mathematica 6.0, I get the same result regardless of whether or how I specify the domain of x.
This is why research should not be done with proprietary mathematical or statistical software any more than lab work should be done by leprechauns.
Off-topic, but I thought it deserved some mention on your blog:
Google and the Wisdom of Clouds
(Business Week Cover Story, 12/13/2007)
Google class debuts at the UW
(Seattle P-I, 2/20/2007)
I have trouble seeing what the problem is. Isn’t Sqrt[x] defined to give a number whose square is x? I mean, it would be different Mathematica was claiming that Sqrt[x^2]=x, but it’s not:
Simplify[Sqrt[x]^2] –> x
Simplify[Sqrt[x^2]] –> Sqrt[x^2]
Simplify[Sqrt[x^2], x ∈ Reals] –> Abs[x]
All of these statements seem correct to me, unless I missed something about numbers having a non-unique square.
I’m going to go out on a limb and argue that Mathematica is right and give two reasons.
Firstly, we are talking about a function of x and requiring that x belongs to the reals. That is, the domain of f(x) is the reals. This does not mean that f(x) can not take values that are not real. In math-speak: the function is a mapping f:x->y from a set x (the domain) to a set y (the range), but there is not reason for these to be the same. So, there is no problem with sqrt(x) taking complex values even if we restrict x to real values. Hence, sqrt(x)^{2}=x because there is no problem with sqrt(x)=ix for negative x.
Secondly, even requiring that the range=domain there is still no problem. To create the supposed bug, we need some peculiar constraint! This is because the problem with complex number crops up at an intermediate stage of the calculation. We have something like f(x)=h(g(x)), where the domain and range of f is. The only objections we can make to f giving negative values are:
i) that g must have a real range; or
ii) that h must have a real domain.
But neither of these are restrictions on f(x). Also, I think the only way to put them into Mathematica would be to do in by hand in some really ad hoc fashion!
IMHO, that is “forgivable” Mathematica behavior, because any other behavior would make Simplify[] almost useless.
However, Mathematica does sometimes commit “unforgivable” errors. For example, floating-point evaluation of SingularValueDecomposition[] sometimes fails with an error message (which is undesirable behavior), and sometimes returns an answer that is grossly wrong (which is unacceptable).
This annoying behavior has persisted for several years. Wolfram Research knows of it, but blames bugs in LAPACK that Wolfram researchers have been unable to diagnose and fix.
However, perhaps this bug has been fixed in version 6 … maybe someone can tell me?
http://faculty.washington.edu/sidles/mma_svd_failure_example/
I use MAXIMA, and whenever it does something that is wrong I figure it’s my punishment for using a free tool instead of shelling out beaucoup bucks for Mathematica.
I end up asking it to do two things and it fails miserably at them. One is to solve say, 6 simultaneous quadratic equations in 6 complex unknowns. You’d think it could do this, but nooooooo.
The other is factorization of nasty equations from GR to simplify them. I ended up having to do most of the work on my own, and then use MAXIMA to multiply it back together and see if it gets the same thing. As it turns out, when you make a sign error or get a factor of two wrong in GR, you end up with orbital simulations that do stuff like violate conservation of energy or angular momentum so you really don’t need symbolic manipulation to tell you that it’s wrong.
If math were a democracy, I’d vote for Eric. Course, Eric would clearly lose to popular opinion.
(And yes, mathematic is perfectly correct here.)
You’d think [Maxima] could do this, but nooooooo.
Nooooooo … because it’s NP-complete!
A nice textbook that explains why is Cox, Little, and O’Shea, Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra.
I haven’t used this particular piece of software, but I’d assume it would also output “1″ for “simplify[x/x, x \in Reals]“, and that no one’s complaining about that.
The consensus is that Mathematica’s simplification Sqrt[x]^2 -> x is correct, in the sense of being “more easy to understand and apply correctly than any reasonable alternative convention”.
But Dave’s question opens the door to a class of deeper questions, that have great importance for our new century, basically involving who owns—or should own—ideas, equations, artworks, and software?
Which is closely linked to the questions, who takes responsibility for teaching, curating, and applying these ideas … and fixing them when mistakes are found?
Needless to say, this question doesn’t have a simple answer (IMHO anyway) … we are all of us engaged in the complex business of community-building, peace-making, and resource-creating.
This too is characteristic of our twenty-first century. We have more ideas at-hand than previous centuries, but on the other hand, our century’s challenges are greater and more urgent.
Hopefully we will do at least as good a job at community-building, peace-making, and resource-creating as previous centuries! 
Or, the related deeper questions, which Stephen Wolfram asked me at ICCS-2007 in Boston, who owns—or should own—ideas, equations, artworks, and software that were “mined” computationally from the space of all possible ideas [the "ideocosm" of Fritz Zwicky], the space of all possible equations [as with brute force searches of axiom systems and proofs, or my being in 1976 the first person to use the Genetic Algorithm to evolve a solution to an unsolved equation in the science literature], or search through spaces of algorithms (as patented by Koza a few years after my PhD research].
Mozart claimed to have plucked music effortlessly from the air, at times, without feeling that he’d written it. Ramanujan claimed to have the goddess Nagiri give him equations. Who owns Mozart’s music, or Ramanujan’s equations, some of which we’ve only understood in the past decade or so?
I think Mathematica’s correct on this.
Sqrt[x], like almost all Mathematica functions, is defined as a single-valued, complex-valued function of a complex variable. It has a well-defined, and documented, choice of branch cut: the -ve real axis, with Sqrt[-1]=i (i.e. the square roots of all negative reals having arguments of pi/2). Since evaluating Sqrt[x]^2 with any x whatsoever will always give x, it’s correct to simplify the expression to just x. That x is restricted to some subset of the complex numbers, such as the reals, doesn’t change that.
Thanks John Sidles, cool to know. Buy hey, even though solving simultaneous quadratic equations may be NP-complete, the problem of solving 6 of them can certainly be done in finite time. I do them regularly, LOL. My complaint is how much time they take. The MAXIMA software won’t even solve 3 of them. Example problem, I, J, K complex:
J = K K + 2 I J,
K = J J + 2 I K,
I = I I + 2 J K.
One of the cool things about systems of simultaneous quadratic equations is that they can do things that are far more bizarre than systems of simultaneous linear equations. These arise in analysis of non perturbative QFT bound states.
Can Mathematica solve the above? There are exactly 8 solutions; they are given by the 8 circulant idempotent 3×3 complex matrices. If someone wants to get started, I’ll contribute one of them: I=J=K=0.
Mathematica computes the Groebner basis as {-k + 729*k^7, j – 81*k^5, -k + 2*i*k + 9*k^4, -i + i^2 + 162*k^6}, from which eight closed-form solutions can be read-off … two solutions with k=0, and one solution for each of the six roots of 729*k^6 = 1.
Pretty much any symbolic manipulation package will compute the Groebner basis of such a small set of equations. By the way, I definitely am *not* an expert in Groebner bases!
In[1]:= Solve[{J == K K + 2 A J, K == J J + 2 A K, A == A A + 2 J K}, {J, K, A}]//InputForm
Out[1]:= {{J -> -1/3, A -> 2/3, K -> -1/3}, {J -> 0, A -> 0, K -> 0}, {J -> 0, A -> 1, K -> 0}, {J -> 1/3, A -> 1/3, K -> 1/3}, {J -> -(-1)^(1/3)/3, A -> 1/3, K -> (-1)^(2/3)/3}, {J -> (-1)^(1/3)/3, A -> 2/3, K -> -(-1)^(2/3)/3}, {J -> -(-1)^(2/3)/3, A -> 2/3, K -> (-1)^(1/3)/3}, {J -> (-1)^(2/3)/3, A -> 1/3, K -> -(-1)^(1/3)/3}}
The point of the above being, that Mathematica (and pretty any other symbol-manipulating program) can find simple solutions to simultaneous algebraic equations if and only the associated Groebner basis is simple enough for its roots to be found “by hand”.
When this happens, there is usually some underlying geometric, algebraic, or informatic symmetry. Often-times the key to the problem is identifying and exploiting this symmetry.
Conversely, if you are solving a system of algebraic problem that has no underlying symmetries whatsoever, then you are tackling a problem that is generically NP-complete … and your chances of exhibiting a simple solution are very poor.
In such cases it is a good idea to ask “Why do I think this problem is important?”
This question *does* have good answers, e.g., “Because I am studying and/or simulating an evolved system.” Examples are DNA sequences, rings of Saturn, immunological evolution, mutual fund prices, stellar clusters, climatology, neurophysiology, turbulent fluid dynamics, etc. — any system that can exhibit behavior of essentially unlimited informatic complexity.
Mathematical disciplines like game theory, control theory, error correction (both classical and quantum), and image reconstruction fall (loosely) into this class too — in the sense that many natural problems in these disciplines are NP-complete.
That’s my 2¢, anyway.
I came across one in using Matlab, which still exists in the latest edition.
Try this:
floor(1.9/0.1)
and this:
floor(3.8/0.1)
You know, float is something very very troublesome.
If Mathematica’s simplify function is anything like its grand-daddy Macsyma’s simplify function there are probably dozens, maybe even hundreds, of flags controlling the meaning of the verb simplify. You probably need to iterate over all the subsets of the set of relevant flags until you get the answer you want. Some things never change.
Dave, I still think Mathematica is correct when it gives you:
Simplify[Sqrt[x]^2, Sqrt[x] in Reals] -> x
and indeed
Simplify[Sqrt[x]^2] -> x
Since the second simplification is true (for all x), the first one must be as well; adding an assumption can only ever enable extra simplifications, it can’t prevent any unconditional ones from applying! All that assumptions ever do is enable some conditional transformations; they don’t redefine what any functions actually mean.
I think what you’re complaining about is that Simplify[Sqrt[x]^2, Sqrt[x] in Reals] is returning an unconditional expression “x” rather than imposing the assumption “Sqrt[x] in Reals” on x itself, to require x ≥ 0. But that’s not how assumptions work, or are advertised as working; the assumption isn’t meant to be incorporated permanently into the simplified expression.
Note that you do get:
Simplify[Sqrt[x] in Reals, x ≥ 0] -> True
but I must confess to being disappointed that you don’t get:
Simplify[x ≥ 0, Sqrt[x] in Reals] -> True
December 29, 2007
December 29, 2007
December 29, 2007
December 29, 2007
December 29, 2007
December 29, 2007
December 29, 2007
December 29, 2007
December 30, 2007
December 30, 2007
December 30, 2007
December 31, 2007
December 31, 2007
December 31, 2007
December 31, 2007