bugGNU Octave - Bugs: bug #40668, -NA is inconsistent with other...


bug #40668: -NA is inconsistent with other operations on NA

Submitter:  Daniel Kraft <domob>
Submitted:  Fri 22 Nov 2013 07:39:27 AM UTC
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Wed 16 Aug 2017 01:27:13 PM UTC, comment #15: 

One of the advantages to switching to checking only the low word for the magic value 1954 is that the sign bit, which is in the upper word, will be ignored.  Thus, -NA and +NA would be the some thing to the isna function which would resolve this bug report in addition to requiring one less equality comparison operator.

Rik <rik5>
Group administrator
Wed 16 Aug 2017 01:23:20 PM UTC, comment #14: 

Study might be too strong a word, but it's definitely possible to use R from C or C++.  The performance details, whether it can be a straight library or whether it is an embedded interpreter are a little vague without doing further research.

Rik <rik5>
Group administrator
Wed 16 Aug 2017 10:48:28 AM UTC, comment #13: 


Good that you have studied this. It is true that R is
THE dominant statistics system. Is it feasible to get
access to their library functions from Octave? If so,
this would be very useful.

Michael Godfrey <godfrey>
Group Member
Wed 16 Aug 2017 05:46:29 AM UTC, comment #12: 

I found the section in the R manual that you did with a description of NA.

The representation of the special values for R numeric and complex types is machine-dependent, and possibly also compiler-dependent. The simplest way to make use of them is to link an external application against the standalone Rmath library which exports double constants NA_REAL, R_PosInf and R_NegInf, and include the header Rmath.h which defines the macros ISNAN and R_FINITE.

If that is not possible, on all current platforms IEC 60559 (aka IEEE 754) arithmetic is used, so standard C facilities can be used to test for or set Inf, -Inf and NaN values. On such platforms NA is represented by the NaN value with low-word 0x7a2 (1954 in decimal).

In liboctave/util/lo-ieee.cc there is both a new and an old NA routine.

__lo_ieee_is_NA (double x)
  lo_ieee_double t;
  t.value = x;
  return (__lo_ieee_isnan (x) && t.word[lo_ieee_hw] == LO_IEEE_NA_HW
          && t.word[lo_ieee_lw] == LO_IEEE_NA_LW) ? 1 : 0;

__lo_ieee_is_old_NA (double x)
  lo_ieee_double t;
  t.value = x;
  return (__lo_ieee_isnan (x) && t.word[lo_ieee_lw] == LO_IEEE_NA_LW_OLD
          && t.word[lo_ieee_hw] == LO_IEEE_NA_HW_OLD) ? 1 : 0;

And finally from lo-ieee.h

#define LO_IEEE_NA_LW_OLD 1954
#define LO_IEEE_NA_HW 0x7FF840F4
#define LO_IEEE_NA_LW 0x40000000
#define LO_IEEE_NA_FLOAT   0x7FC207A2
#define LO_IEEE_NA_HW_OLD 0x7ff00000

The old and new are weird, and from the R documentation I think Octave is trying to be too specific.  The code shouldn't need to test the upper word.  I think testing only for isnan and the lower word equal to 1954 is enough.  But this brings up another issue.  Like the code which determines endianness, this code writes to one struct member, and then reads from another, the behavior of which is not guaranteed in C++.  Either this file needs to be a straight C file and compiled with a C compiler, or it should be re-coded.

Here's an attempt

__lo_ieee_is_NA (double x)
  return ((__lo_ieee_isnan (x) && (reinterpret_cast<uint32_t *> (&x) == 1954)) ? 1 : 0;)

I think this depends on a little-endian machine though.

Rik <rik5>
Group administrator
Tue 15 Aug 2017 11:36:58 PM UTC, comment #11: 

The R manual I downloaded has

NA is a logical constant of length 1 which contains a missing value indicator.
NA can be coerced to any other vector type except raw. There are also constants
NA_integer_, NA_real_, NA_complex_ and NA_character_ of the other atomic vector
types which support missing values: all of these are reserved words in the R

I searched for 'is represented by', but couldn't find anything.  My motivation was partly to pick off an easy bug, although I like the idea of NA as being more explicit than NaN about Not Available.

Matlab isn't particularly known for their statistics package, and likewise Octave.  If we wanted too we could outsource statistics functions to their library and thereby get some very high quality code.  In that case, it might be nice to support NA.

Rik <rik5>
Group administrator
Tue 15 Aug 2017 06:15:40 PM UTC, comment #10: 

The R manual says:

On such platforms @code{NA} is represented by
the @code{NaN} value with low-word @code{0x7a2}
(1954 in decimal).

The file src/main/arithmetic.c has this:

static double R_ValueOfNA(void)
    /* The gcc shipping with Fedora 9 gets this wrong without
     * the volatile declaration. Thanks to Marc Schwartz. */
    volatile ieee_double x;
    x.word[hw] = 0x7ff00000;
    x.word[lw] = 1954;
    return x.value;

I don't think R has single precision floating point values, so there probably isn't a single-precision NA.  I don't see one, anyway.

Originally we used that value, but it looks like it was changed, possibly to accomodate single-precision values.

At this point, I'd be glad to drop NA values.  In more than 15 years, we've never done anything significant with them.  It's not like all the stats functions in Octave actually pay attention to NA.

John W. Eaton <jwe>
Group administrator
Tue 15 Aug 2017 05:31:45 PM UTC, comment #9: 

This is a very small bug, but it does seem reasonable to expand __lo_ieee_is_NA to recognize both NA and -NA.

Is there a good way to get the value of NA and -NA from R?

Mike did a test back in 2016.  Maybe it is enough to write out the binary values and then look at the file in a binary editor.  If I had a copy of the file I would be willing to take a look.

Rik <rik5>
Group administrator
Wed 16 Nov 2016 11:17:02 PM UTC, comment #8: 

I assume you are asking about unpatched Octave as we know it? That's what I tested here.

I know almost nothing about R, but here is what I tried:

$ R
> fid = file ("foo.dat", "wb")
> writeBin (c (1, -1) * NA, fid)
> close (fid)

$ octave
>> fid = fopen ("foo.dat", "rb");
>> x = fread (fid, "double");
>> fclose (fid);
>> x
x =


So that's one direction. In the other direction:

$ octave
>> fid = fopen ("foo.dat", "wb");
>> fwrite (fid, [1, -1] * NA, "double");
>> fclose (fid);

$ R
> fid = file ("foo.dat", "rb")
> x = readBin (fid, double (), n = 2)
> close (fid)
> x
[1] NaN NaN

Is that what you are looking for? If not what can I test differently?

Mike Miller <mtmiller>
Group Member
Mon 25 Nov 2013 03:30:12 AM UTC, comment #7: 

Could someone try the following and see if it works?

Write NA and -NA to a file in binary format in R and read the file in Octave.  Do you get NA for both?

Write NA and -NA to a file in binary format in Octave and read the file in R.  Do you get NA for both?

Being able do do this (at least for NA, not sure about -NA) was the original motivation for having NA in Octave implemented like the NA from R.  So if it is not possible to do this anymore, we also need to find a way to restore this capability.

John W. Eaton <jwe>
Group administrator
Mon 25 Nov 2013 02:34:21 AM UTC, comment #6: 

> Are there possible problems that can be introduced by modifying __lo_ieee_is_NA() accordingly?

Well, the point was to do whatever R was doing. What does R do internally with -NA? In its REPL, -NA is still NA.

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Fri 22 Nov 2013 03:36:47 PM UTC, comment #5: 

As I understand it, doesn't IEEE 754 allow for distinct (quiet) NaN values, of which NA is a particular one?  Those values should be propagated through all computations, so as long as you don't perform operations between NA and NaN but only between NA and ordinary numbers, NA should be preserved (also for external codes that "just" do calculations with the numbers passed in).

The sign bit is a bit different from the NaN payload bits, and thus seemingly the unary minus flips the sign bit even though the payload itself stays the same as it should.  Based on this interpretation, it makes sense to me to ignore the sign bit when checking for NA.  However, I understand that it is probably best not to rely on this behaviour in one's code (and I have changed mine already accordingly).

Are there possible problems that can be introduced by modifying __lo_ieee_is_NA() accordingly?

Daniel Kraft <domob>
Fri 22 Nov 2013 02:09:27 PM UTC, comment #4: 

The particular values for NA were taken from R and were chosen so that you could exchange binary data files between Octave and R and preserve NA.

But in more than a decade, nothing else about NA has really been implemented in Octave.  And I'm not sure we should start now.  It's pretty much impossible for us to guarantee that all operations on NA return NA rather than NaN.  It's certainly impossible for us to guarantee that for operations that involve external numerical libraries like ATLAS or OpenBLAS.

John W. Eaton <jwe>
Group administrator
Fri 22 Nov 2013 01:16:17 PM UTC, comment #3: 

I wondered about the performance hit of doing it higher up as well. Does your quick hack at least fix the problem?

I also don't know how the specific NA value was chosen and what the implications might be for considering it equal ignoring the sign bit. Was it selected because it's a value that floating point math operations will never evaluate to? Is the same true for the negative value?

Mike Miller <mtmiller>
Group Member
Fri 22 Nov 2013 08:44:06 AM UTC, comment #2: 

The problem here seemingly is that the unary minus changes the sign bit, which makes Octave no longer recognise its special NaN-pattern for NA.  A possible fix is to change __lo_ieee_is_NA to check for both patterns (the original one and the one with inverted sign bit), like this:

__lo_ieee_is_NA (double x)
#if defined (HAVE_ISNAN)
  lo_ieee_double t;
  t.value = x;
  return (isnan (x) && (t.word[lo_ieee_hw] == LO_IEEE_NA_HW || t.word[lo_ieee_hw] == 0xFFF840F4)
          && t.word[lo_ieee_lw] == LO_IEEE_NA_LW) ? 1 : 0;
  return 0;

(This is of course not a properly done patch, just a quick hack to try it out.)

Since I'm not an expert with IEEE floating-point numbers and things like that:  Is this the correct way to fix this bug, or should we rather add a check to unary minus operations so that they don't invert the sign bit in case of NA?  This, however, seems to me as if it would hurt performance for each and every unary minus operation, and would also possibly have to be done in all the different implementations of unary minus (for scalars, complex values, matrices, sparse matrices and so on).  Thus to me the change of isna() seems better.

If this sounds like a good approach, I will write a nicer patch (including some comments) and provide a hg commit diff.

Daniel Kraft <domob>
Fri 22 Nov 2013 07:40:35 AM UTC, comment #1: 

I will try to work on a patch.

Daniel Kraft <domob>
Fri 22 Nov 2013 07:39:27 AM UTC, original submission:  

Currently, -NA evaluates to NaN, while a lot of other operations on NA result in NA.  In particular:

+NA, 0 - NA, -1 NA, Inf NA, 0 * NA, NA / 0, 0 / NA, and others (especially of course less edge-case like operations, too).

See also https://mailman.cae.wisc.edu/pipermail/help-octave/2013-November/061056.html.  I propose that -NA should evaluate to NA to make the behaviour more consistent.

Daniel Kraft <domob>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by domob (Submitted the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.


    Follow 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-11-16 mtmiller Release3.6.4 dev
    2013-11-22 mtmiller Item GroupOther Incorrect Result
        StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-ee12.
    Corresponding source code