9
$\begingroup$

I'm working in Mathematica 13.0 and I type

a=0.00001;
IntegerPart[a/a]

and well, the output is 0.

Try with almost any other value of a, say a=0.01, and the result is correctly 1. I did some trials and I found that it fails again with a=0.000000001. I also found that a possible fix is

a=0.00001 // Rationalize;

(then I got tired). What is going on? Can this have to do with the MachinePrecision?

More in general, how do I fix it once and for all for the entire Notebook?

$\endgroup$
7
  • 4
    $\begingroup$ RealDigits[a/a] => {{9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9}, 0} In the memory this is a number less than 1. Yes it has to do with machine epsilon $\endgroup$ Commented Aug 8, 2022 at 14:31
  • 2
    $\begingroup$ Also interesting example: arr = RandomReal[1, 100000]; Tally[RealDigits /@ (arr/arr)] => {{{{9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9}, 0}, 15338}, {{{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1}, 84662}}. ~15% numbers giving not 1 when you try to divide number on itself $\endgroup$ Commented Aug 8, 2022 at 14:38
  • 4
    $\begingroup$ Try IntegerPart[Divide[a, a]]. It has been noted before on this site that "a / a" is parsed as Times[a, Power[a, -1]]. The roundoff error in 1 / a = 99999.99999999999 means a * (1/a) does not cancel (and sometimes doesn't for other values of a). $\endgroup$
    – Michael E2
    Commented Aug 8, 2022 at 14:42
  • 4
    $\begingroup$ @KirillBelov You might be interested in this: Table[a = 10.`16^-k; IntegerPart[a/a], {k, 100}]. The use of arbitrary-precision means more bits than in machine precision are used to represent the floating-point numbers (plus precision tracking, which is unimportant here). All the extra precision does is change where the problem occurs. $\endgroup$
    – Michael E2
    Commented Aug 8, 2022 at 14:58
  • 2
    $\begingroup$ See also this and this. $\endgroup$
    – user293787
    Commented Aug 9, 2022 at 4:31

1 Answer 1

11
$\begingroup$

[H]ow do I fix it once and for all for the entire Notebook?

I suppose it depends on how you want to treat results that contain round-off error. You cannot really get rid of the problem, only shift where the problem arises. For instance, in IntegerPart[x], IntegerPart (nor usually the program/programmer) knows at that point whether the roundoff error that led to x is positive, negative, or zero. If you fix IntegerPart[x] when x comes from a/a and is one bit less than 1., then it will be broken when another x has a roundoff error that increases it to one bit less than 1.. One may be willing to live with it being broken in this way. Here's a way based on Equal comparing with tolerance:

intPart // ClearAll;
intPart[x_?NumericQ] /; x == Round[x] := Round[x];
intPart[x_?NumericQ] := IntegerPart[x];

intPart[a/a]
(*  1  *)

Make it Listable if your code relies on IntegerPart being listable. You can change the tolerance used by Equal with Internal`$EqualTolerance.

Another approach is to "correct" real numbers that are close to integers:

snapToInteger // ClearAll;
snapToInteger[expr_] := 
 expr /. {x_Real /; x == Round[x] :> 
    SetPrecision[Round[x], Precision[x]]};

IntegerPart[snapToInteger[a/a]]
(*  1  *)

Either approach is likely to incur a performance hit, which may or may not be important.

For what it's worth, this works on x = a/a, because of the nature of the rounding error and how MatchQ works on machine reals and bignums, and it's a bit faster:

intPart // ClearAll;
intPart[x_Real] /; MatchQ[x, N[Round[x], Precision[x]] := 
  Round[x];
intPart[x_?NumericQ] := IntegerPart[x];
$\endgroup$
2
  • 1
    $\begingroup$ I have to say that using Divide[a/a], as you suggested under my question, is also (often) a quick and valid fix. $\endgroup$
    – sonarventu
    Commented Aug 10, 2022 at 9:10
  • 1
    $\begingroup$ @DavideVenturelli You may have noticed that if x is a symbol (without a numeric or other value), Divide[x, x] evaluates to x/x, that is, Times[x, Power[x, -1]]. Sometimes it's not so easy to replace / by Divide[], and since there was no response to my comment, I thought it might not work in all your use-cases. $\endgroup$
    – Michael E2
    Commented Aug 10, 2022 at 13:18

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.