4
$\begingroup$

Consider this code

x = 0.109354682484;
IntegerPart[x/(x/2)]
(* 1 *)

Precision[x]
(* MachinePrecision *)

Why does it give 1 ?

Version number: 9.0 on Mac 10.9.2

screenshot

enter image description here

Update:

If we use an undefined variable, IntegerPart[x0/(x0/2)] gives 2. Since Mathematica never gives warnings about this x0, I'm assuming for any x0 it is true.

If we calculate the same integer part using fortran, we get 2 instead of 1.

program main
implicit none

real(8):: x=0.109354682484
real(8):: y=1.4

write(*,*) int(x/(x/2))
write(*,*) nint(x/(x/2))
write(*,*) int(y)

end program main

compiled with ifort -O0 main.f90

output of above fortran code is

       2
       2
       1

according to here, int is a fortran intrinsic function that calculate the integer part.

Is this a bug?

Improve this question
$\endgroup$
6
  • 3
    $\begingroup$ Look at the FullForm of the result of x/(x/2). $\endgroup$
    – chuy
    Commented Apr 1, 2014 at 17:23
  • $\begingroup$ I'm aware of that, but why fortran gives 2, with the same precision? $\endgroup$
    – xslittlegrass
    Commented Apr 1, 2014 at 17:46
  • $\begingroup$ Then I guess it comes down to what you want IntegerPart[1.9999999999999998] to return? $\endgroup$
    – chuy
    Commented Apr 1, 2014 at 17:50
  • 3
    $\begingroup$ Does your Fortran compiler optimize the expression x/(x/2)? $\endgroup$
    – Michael E2
    Commented Apr 1, 2014 at 18:00
  • $\begingroup$ @MichaelE2 I only use -O0 option but I don't know whether there are other optimizations. Are you aware of other optimizations? $\endgroup$
    – xslittlegrass
    Commented Apr 1, 2014 at 18:07

1 Answer 1

Reset to default
7
$\begingroup$

The explanation is interesting here. I tried the same in C++, and worked a bit extra to make sure the compiler won't optimize away the divisions (looking at the assembly output, it may optimize it away if you're not careful). Indeed, I get 2 with C++.

And here's why:

C++ does the equivalent of 

IntegerPart@Divide[x, Divide[x, 2]]

(* ==> 2 *)

while in Mathematica you're computing the equivalent of

x*(1/((1/2)*x))

which gives a different result.

If I do the calculation as (x* (1/ (0.5*x))) in C++, then the result is less than 2 and the integer part is 1.

Relevant reading:

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks! I never knew that a/b is interpreted as a*(1/b). That makes me worried a little bit about all the codes I have written. Actually I caught on this because of a strange behavior of my package, and it took me a whole morning to track down to problem to IntegerPart I used in that package. $\endgroup$
    – xslittlegrass
    Commented Apr 1, 2014 at 18:54
  • 1
    $\begingroup$ The observed behavior is wrong under IEEE-754 as well as generic binary floating-point implementations. The result of x/2 or 0.5 * 2 has the same mantissa as x, and the exponent is decreased by one. Dividing two numbers having the same mantissa then results in a mantissa of 1 (stored implicitly), and the exponent is the difference in operand exponents. The only way for a binary implementation to not get exactly 2 for the result of x/(x/2) or x/(0.5*x) is if the scaled value becomes denormal. $\endgroup$
    – Ben Voigt
    Commented Apr 1, 2014 at 23:13
  • $\begingroup$ @BenVoigt But as I said it does not compute x/(0.5*x). It computes x * (1/(0.5*x)). I don't know what's correct by the standard, but the result is precisely the same as what I get from C++, which just uses a few simple floating point instructions (based on the asm output). Mathematica also uses hardware floating point arithmetic by default. $\endgroup$
    – Szabolcs
    Commented Apr 2, 2014 at 4:55
  • $\begingroup$ @BenVoigt The reason why it does the arithmetic in this order is explained in the post I linked to. $\endgroup$
    – Szabolcs
    Commented Apr 2, 2014 at 4:57
  • $\begingroup$ Oh right, I thought about only one reciprocal and there are two. $\endgroup$
    – Ben Voigt
    Commented Apr 2, 2014 at 6:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.