Why does NaN^0 == 1

This is referenced in the help page referenced by ?'NaN'

"The IEC 60559 standard, also known as the ANSI/IEEE 754 Floating-Point Standard.

http://en.wikipedia.org/wiki/NaN."

And there you find this statement regarding what should create a NaN:

 "There are three kinds of operations that can return NaN:[5]
       Operations with a NaN as at least one operand.

It is probably is from the particular C compiler, as signified by the Note you referenced. This is what the GNU C documentation says:

http://www.gnu.org/software/libc/manual/html_node/Infinity-and-NaN.html

" NaN, on the other hand, infects any calculation that involves it. Unless the calculation would produce the same result no matter what real value replaced NaN, the result is NaN."

So it seems that the GNU-C people have a different standard in mind when writing their code. And the 2008 version of ANSI/IEEE 754 Floating-Point Standard is reported to make that suggestion:

http://en.wikipedia.org/wiki/NaN#Function_definition

The published standard is not free. So if you are have access rights or money you can look here:

http://ieeexplore.ieee.org/xpl/mostRecentIssue.jsp?punumber=4610933

The answer can be summed up by "for historical reasons".

It seems that IEEE 754 introduced two different power functions - pow and powr, with the latter preserving NaN's in the OP case and also returning NaN for Inf^0, 0^0, 1^Inf, but eventually the latter was dropped as explained briefly here.

Conceptually, I'm in the NaN preserving camp, because I'm coming at the issue from viewpoint of limits, but from convenience point of view I expect current conventions are slightly easier to deal with, even if they don't make a lot of sense in some cases (e.g. sqrt(-1)^0 being equal to 1 while all operations are on real numbers makes little sense if any).

Yes, I'm late here, but as R Core member who was involved in this design, let me recall what I commented above. NaN preserving and NA preserving work "equivalently" in R, so if you agree that NA^0 should give 1, NaN^0 |-> 1 is a consequence.

Indeed (as others said) you should really read R's help pages and not C or IEEE standards, to answer such questions, and SimonO101 correctly cited

1 ^ y and y ^ 0 are 1, always

and I'm pretty sure that I was heavily involved (if not the author) of that. Note that it is good, not bad, to be able to provide non-NaN answers, also in cases other programming languages do differently. The consequence of such a rule is that more things work automatically correctly; in the other case, the R programmer would have been urged to do more special casing herself.

Or put differently, a simple rule as the above (returning non-NaN in all cases) is a good rule, because it propagates continuity in a mathematical sense: lim_x f(x) = f(lim x). We have had a few cases where it was clearly adavantageous (i.e. did not need special casing, I'm repeating..) to adhere to the above "= 1" rule, rather than to propagate NaN. As I said further up, the sqrt(-1)^0 is also such an example, as 1 is the correct result as soon as you extend to the complex plane.