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