Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.math.mathematica

Topic: Mathematica not IEEE-754 compliant?
Replies: 4   Last Post: Nov 6, 2013 12:21 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Alain Cochard

Posts: 7
Registered: 5/6/13
Mathematica not IEEE-754 compliant?
Posted: Oct 28, 2013 11:14 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply



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@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




Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.