From: Ted Harding <Ted.Harding_at_manchester.ac.uk>

Date: Sat 20 Jan 2007 - 19:53:07 GMT

*> [2,] 0.7071068
*

*> [3,] 0.0000000
*

> [1] 1

> [1] FALSE

*>
*

*> I realize that it has something to do with the conversion of the square
*

*> roots into decimals. But shouldn't it then give me some number other
*

*> than 1 as the result for sum(A3^2)? Are there other ways to do this
*

*> than what I've tried? I'm trying to confirm that A3 is a unit vector.
*

E-Mail: (Ted Harding) <Ted.Harding@manchester.ac.uk> Fax-to-email: +44 (0)870 094 0861

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Jan 21 06:59:04 2007

Date: Sat 20 Jan 2007 - 19:53:07 GMT

On 20-Jan-07 Robert Barber wrote:

*> Hi,
**>
*

> I apologize if there is a simple answer to this question that I've

*> missed. I did search the mailing list but I might not have used the
**> right keywords. Why does sum(A3^2) give the result of 1, but
**> sum(A3^2)==1 give the result of FALSE?
**>
*

>> A3<-matrix(nrow=3,c(1/(2^.5),1/(2^.5),0)) >> A3

> [,1]

> [1,] 0.7071068

>> sum(A3^2)

> [1] 1

>> sum(A3^2)^.5

> [1] 1

>> sum(A3^2)==1 # here's the part I don't understand

> [1] FALSE

>> sum(A3^2)^.5==1 # here's the part I don't understand

> [1] FALSE

This is an instance of what must be a candidate for the MFAQAT (most frequently asked qustion of all time).

The nub of the matter can be found in FAQ 7.31:

http://www.r-project.org/ --> FAQs

where, at 7.31, it says:

7.31 Why doesn't R think these numbers are equal?

The only numbers that can be represented exactly in R's numeric type are integers and fractions whose denominator is a power of 2. Other numbers have to be rounded to (typically) 53 binary digits accuracy. As a result, two floating point numbers will not reliably be equal unless they have been computed by the same algorithm, and not always even then. For example

R> a <- sqrt(2) R> a * a == 2 [1] FALSE R> a * a - 2 [1] 4.440892e-16

The function all.equal() compares two objects using a numeric tolerance of .Machine$double.eps ^ 0.5. If you want much greater accuracy than this you will need to consider error propagation carefully.

For more information, see e.g. David Goldberg (1991),
"What Every Computer Scientist Should Know About Floating-Point
Arithmetic", ACM Computing Surveys, 23/1, 5-48, also
available via

http://docs.sun.com/source/806-3568/ncg_goldberg.html.

However, what this FAQ does not point out is that "==" tests for exact equality, and even the 'help' page

?"=="

is not as explicit about this:

For numerical values, remember '==' and '!=' do not allow for the finite representation of fractions, nor for rounding error. Using 'all.equal' with 'identical' is almost always preferable. See the examples.

Nor does the help

?all.equal

give the explanation:

'all.equal(x,y)' is a utility to compare R objects 'x' and 'y' testing "near equality". ...

and a user who is not aware of the imprecision problems inherent in computer arithmetic may well not see the point of testing for "near equality" when the user expects exact equality mathematically.

Statistically speaking, the frequency of occurrence of variants of this question is a phenomenon!

I hypothesise that it arises from a combination of two main circumstances:

- Many users who are unfamiliar with the technical details of finite-length binary arithmetic will not be expecting that there could be a problem of this kind in the first place. So, when it occurs, they will simply be puzzled.
- It's actually quite difficult to find your way to the above explanation in the FAQs. First, you need to anticipate that this is the sort of thing that will be a FAQ. If you're subject to (a) above, you may not be thinking on these lines.

Secondly, even if you get as far as looking at the FAQs at the above URL, you have a lot of scrolling down to do before you find the question

7.31 Why doesn't R think these numbers are equal?

and even then, if you blink you'll miss it: trying it just now, I in fact didn't spot it in the list of questions immediately below the header

7 R Miscellanea

even though I already knew it was there somewhere and was actively looking for it. It was only when I landed on the full FAQ itself that I recognised it.

It would have been quicker (I just tried that too) to start at the top of the FAQs page, and use the browser's Search tool to search for == : the sixth occurrence of "==" was it!

But, even then, you still need to be thinking that the answer is to be found in connection with "==". If you're subject to (a), you'll be wondering instead why the answer was wrong.

This is such a frequently occurring issue that I feel there is a case for a prominently displayed link, maybe in FAQs but maybe better at the top level of r-project.org, to a brief but adequate discourse with a title like

Arithmetic [im]precision and apparent errors in R

What do others think?

Ted.

E-Mail: (Ted Harding) <Ted.Harding@manchester.ac.uk> Fax-to-email: +44 (0)870 094 0861

Date: 20-Jan-07 Time: 19:53:03 ------------------------------ XFMail ------------------------------ ______________________________________________R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Jan 21 06:59:04 2007

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sun 21 Jan 2007 - 02:30:22 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*