2
$\begingroup$

The compile function I have written computes a table of relative frequencies of some data. Thus, the output is a list of increase elements, p, in [0,1]. Since I have say 1000 data points the elements p are equal to x/1000 and x in [0,1000]. Hence, each p should have a maximum of three positions after the decimal point. However, if I run the compile function then values having many more positions after the decimal point appear. However, this appears to be true not for all elements of the resulting list.

For example, let the true p equal 0.894, then the compile output reads p = 0.8949999999999.

How can I force the result to be 0.894?

Unfortunately, Round[p,1./1000.]does not work within the compile function. It does not work either if I apply Roundto the resulting list. What works, but seem inefficient, is Round[list,1/1000]//N. Thanks for help.

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ Executing RealDigits[894/1000, 2] will help you see what's going on: 0.894 has a nontrivial infinitely repeating expansion in base 2 (exactly like, say, 1/3 or 1/7 have repeating expansions base 10), proving that it's impossible to represent 0.894 exactly in base 2 floating point representations. $\endgroup$
    – whuber
    Commented Oct 14, 2012 at 15:23

1 Answer 1

Reset to default
4
$\begingroup$

This is a question about floating-point arithmetic, not Mathematica as such.

The Mathematica part of the answer is that using Compile means that all computations happen with machine numbers, and then it's a question of understanding what floating-point arithmetic does. See also this.

A simple demonstration of the dangers of using floating-point numbers as if they had infinite accuracy:

y = 0.01 $MachineEpsilon
x = 1. + y;
x - 1

Mathematica graphics

Improve this answer
$\endgroup$
8
  • $\begingroup$ I don't like this example. It should give False under any meaningful definition of epsilon, but doesn't because of the extra tolerance of Log[2, 10^Internal`$EqualTolerance] bits applied by Equal. If you want a number for which 1 + x == x gives True and don't wish to unduly abuse IEEE754, I would suggest using x = (1 + 10^Internal`$EqualTolerance)/2 $MachineEpsilon instead. Alternatively, you could use Internal`CompareNumeric[0, 1 + $MachineEpsilon/2, 1]. $\endgroup$
    – Oleksandr R.
    Commented Oct 14, 2012 at 22:55
  • $\begingroup$ @OleksandrR. I don't understand what you're saying I'm afraid. Can you unpack this a bit? (also, any thoughts on this?) $\endgroup$
    – acl
    Commented Oct 14, 2012 at 23:01
  • $\begingroup$ $MachineEpsilon is defined as the smallest number for which x + x $MachineEpsilon == x gives False in an ordinary floating-point comparison. In Mathematica, there is an extra tolerance--Log[2, 10^Internal`$EqualTolerance] bits are effectively chopped off each comparand before the comparison is done, so by writing 1 + $MachineEpsilon == 1 you are basically requesting the result of 1. == 1, which entirely defeats the point. Because of the tolerance, you need to use a larger number than $MachineEpsilon for the example to be meaningful. $\endgroup$
    – Oleksandr R.
    Commented Oct 14, 2012 at 23:09
  • $\begingroup$ @OleksandrR. I see your point (I was just experimentally looking for a small enough number so to my shame did not even know what MachineEpsilon is supposed to be). How do you feel about x = 1. + 0.001 $MachineEpsilon;x - 1 as an example, then? $\endgroup$
    – acl
    Commented Oct 14, 2012 at 23:14
  • $\begingroup$ That'd be better, but 0.001 is too small. Use 0.5 instead. :) Sorry, no thoughts on the tridiagonal matrix construction--just stopped by for 10 minutes as a break from thesis-writing. I'll have a proper look after I'm done, 2 weeks from now. $\endgroup$
    – Oleksandr R.
    Commented Oct 14, 2012 at 23:19

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.