Re: Mathematica not IEEE-754 compliant?
- To: mathgroup at smc.vnet.net
- Subject: [mg131961] Re: Mathematica not IEEE-754 compliant?
- From: Alain Cochard <Alain.Cochard at unistra.fr>
- Date: Wed, 6 Nov 2013 00:37:24 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-outx@smc.vnet.net
- Delivered-to: mathgroup-newsendx@smc.vnet.net
- References: <l4n97h$b25$1@smc.vnet.net> <l4sjo2$pkv$1@smc.vnet.net>
danl at wolfram.com writes:
> There is a subtlety having to do with Mathematica semantics rather
> than IEEE arithmetic. n/d is intertpreted, and evaluated, as
> n*(d^(-1)). If you go through the various computations to check on how
> your example will behave, I suspect it will come out fine, with the
> ELP ("error in last place") discrepancy between n*(d^(-1)) and n/d
> resulting from roundoff error in the process of taking the reciprocal
> and then multiplying.
>
> To get the division semantics directly in Mathematica one can instead use
>
> Divide[n,d]
>
> In this case your example will give the result you had expected.
>
> I should mention that once the actual arithmetic is performed in
> Mathematica, when operands are machine doubles then it is done using
> hardware operations that are invoked from C. If the C compiler used in
> the build process is IEEE compliant, then those operations will be as
> well. To the best of my knowledge all compilers we use to build
> Mathematica are IEEE compliant.
That's a subtlety, all right!
Yes, knowing about it, all my results make sense to me:
- d/n is never more that one machine number away from the closest
machine number (more precisely, for d and n random numbers in [1/2, 1],
it gives the closest machine number about 73% of the time and a number
one machine number away from the closest the remaining 27%);
- Divide[d,n] always gives the closest machine number.
I must confess I had not read the documentation for "/"... But the
subtlety does not seem to be documented there. Indeed, in the Divide
page
http://reference.wolfram.com/mathematica/ref/Divide.html
I read that
x/y or Divide[x, y] is equivalent to x y^-1
The Details section does not say more, it seems. I wonder then if there
is a place where one can learn about this kind of subtleties...
Thanks again for your time.
a.
--
EOST (Ã?cole et Observatoire des Sciences de la Terre)
IPG (Institut de Physique du Globe) | alain.cochard at unistra.fr
5 rue René Descartes [bureau 106] | Phone: +33 (0)3 68 85 50 44
F-67084 Strasbourg Cedex, France | Fax: +33 (0)3 68 85 01 25