Mathematica not IEEE-754 compliant?
- To: mathgroup at smc.vnet.net
- Subject: [mg131909] Mathematica not IEEE-754 compliant?
- From: Alain Cochard <Alain.Cochard at unistra.fr>
- Date: Mon, 28 Oct 2013 23:25:16 -0400 (EDT)
- 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
Hello.
My understanding is that, when machine-precision numbers are used, Mathematica
follows the IEEE-754 standard, hence results should be identical to
other IEEE-754 compliant software+hardware combinations when the +, -,
*, /, operations are used.
Well, I am aware of the fact that floating point arithmetics can be
tricky and that systems that do follow the IEEE-754 standard may give
different results, but it seems to me it should not be the case here.
More precisely, my understanding is that, given any one of these
arithmetic operations, the machine number provided as the result should
be the machine number closest to the exact result (for the default
rounding mode).
I believe I show below a case where this is not true in Mathematica.
NB: the tests have been done for Mathematica versions 7.0 and 8.0 on two
different GNU/Linux machines. This is for the division of two numbers,
'num' and 'den' (for 'numerator' and 'denominator').
----
(* Infinite precision versions of the numbers: *)
numInfPrec=2^^11100011010111111010100100110001100111111111111101111*2^-3;
denInfPrec=2^^11100011010111111010100100110001100111111111111110111*2^-3;
(* Transform to machine-precision numbers: *)
num=1.*numInfPrec;
den=1.*denInfPrec;
(* The strings of 0s and 1s above start with a '1' and are 53-digit
long, hence the above numbers correspond to regular normalized numbers
of the 64-bit long IEEE-754 format -- indeed, the outputs of
RealDigits[num,2] and RealDigits[den,2] are consistent with the above
numbers. *)
res=num/den; (* Actual computation under testing *)
bin=RealDigits[res,2] (* Corresponding binary representation *)
This gives:
{{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 0, 1, 1, 0}, 0}
I contend that is this not the machine number closest to the exact
result and that, instead,
1.1111111111111111111111111111111111111111111111110111 * 2^-1 is the
closest number, called below "presumablyClosestInfPrec" for its
infinite precision version; this number differs from the Mathematica number
above by the last binary digit of the mantissa -- it is the next
machine number.
Checking:
---------
(* Infinite precision number corresponding to the result, res: *)
resInfPrec=FromDigits[bin,2];
presumablyClosestInfPrec=
2^^11111111111111111111111111111111111111111111111110111*2^-53;
resExactInfPrec=numInfPrec/denInfPrec;
(* Difference between Mathematica result and exact result: *)
diffMma=N[resInfPrec-resExactInfPrec,3]//InputForm
(* This gives -1.1022302462516`3.*^-16 *)
(* Diff between presumably closest machine number and exact result: *)
diffPresumablyClosest=N[presumablyClosestInfPrec-resExactInfPrec,3]//InputForm
(* This gives 7.9927783736023861873`3.*^-19 *)
This shows that the exact result of the division is between the two
(adjacent) machine numbers and indeed closer to the presumably closest
number.
----
Another way to illustrate the problem:
It just so happens that num=(1/e)-2 and den=(1/e)-1, where e=1. 10^-15.
It is very easy to check that num/den obtained by Mathematica is different
from num/den obtained with another IEEE-754 compliant software (I have
tried with Matlab, Octave, and GNU Fortran):
Mathematica:
e=1. 10^-15;
num=(1/e)-2 ; InputForm[num] (* gives 9.999999999999979*^14 *)
den=(1/e)-1 ; InputForm[den] (* gives 9.999999999999989*^14 *)
res=num/den ; InputForm[res] (* gives 0.9999999999999989 *)
matlab:
e=1e-15;
num=(1/e)-2 % gives 9.999999999999979e+14
den=(1/e)-1 % gives 9.999999999999989e+14
res=num/den % gives 9.999999999999990e-01
num and den are the same for Mathematica and Matlab, while res=num/den differ.
I have also checked that
(*) the 52+1 binary representations of num and den obtained with
Matlab and Octave (with decimal-to-binary converters), and GNU Fortran
(with the b64.64 format), are identical to the binary reresentations
given by Mathematica;
(*) the 52+1 binaray representations of 'res' obtained with Matlab,
Octave, and GNU Fortran, are all identical to the representation of
what I refer above by the "closest machine number".
----
Hopefully, if Mathematica is not faulty here, someone can point out the
mistake I am making.
Thank you.
Alain Cochard
--
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