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: firstname.lastname@example.org
- Delivered-to: email@example.com
- Delivered-to: firstname.lastname@example.org
- Delivered-to: email@example.com
- References: <firstname.lastname@example.org> <email@example.com>
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