
Mathematica not IEEE754 compliant?
Posted:
Oct 28, 2013 11:14 PM


Hello.
My understanding is that, when machineprecision numbers are used, Mathematica follows the IEEE754 standard, hence results should be identical to other IEEE754 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 IEEE754 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 machineprecision numbers: *) num=1.*numInfPrec; den=1.*denInfPrec;
(* The strings of 0s and 1s above start with a '1' and are 53digit long, hence the above numbers correspond to regular normalized numbers of the 64bit long IEEE754 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[resInfPrecresExactInfPrec,3]//InputForm (* This gives 1.1022302462516`3.*^16 *)
(* Diff between presumably closest machine number and exact result: *) diffPresumablyClosest=N[presumablyClosestInfPrecresExactInfPrec,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 IEEE754 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=1e15; num=(1/e)2 % gives 9.999999999999979e+14 den=(1/e)1 % gives 9.999999999999989e+14 res=num/den % gives 9.999999999999990e01
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 decimaltobinary 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@unistra.fr 5 rue RenÃ© Descartes [bureau 106]  Phone: +33 (0)3 68 85 50 44 F67084 Strasbourg Cedex, France  Fax: +33 (0)3 68 85 01 25

