floats: add functions for handling NaN payloads

This handles the simple case of ignoring sign and only dealing with
quiet NaN values since signalling NaN would require that we check that
at least one of the low bits was set and fail somehow if not.

Note that Go treats all NaNs as silent.

floats/floats.go

@@ -490,6 +490,30 @@ func MulTo(dst, s, t []float64) []float64 {
 	return dst
 }
 
+const (
+	nanBits = 0x7ff8000000000000
+	nanMask = 0xfff8000000000000
+)
+
+// NaN returns an IEEE 754 "quiet not-a-number" value with the
+// payload specified in the low 51 bits of payload.
+// The NaN returned by math.NaN has a bit pattern equal to NaN(1).
+func NaN(payload uint64) float64 {
+	payload &^= nanMask
+	return math.Float64frombits(nanBits | payload)
+}
+
+// NaNPayload returns the lowest 51 bits of an IEEE 754 "quiet
+// not-a-number" and true, or zero and false if f is not NaN
+// or is not quiet.
+func NaNPayload(f float64) (payload uint64, ok bool) {
+	b := math.Float64bits(f)
+	if b&nanBits != nanBits {
+		return 0, false
+	}
+	return b &^ nanMask, true
+}
+
 // Nearest returns the index of the element in s
 // whose value is nearest to v.  If several such
 // elements exist, the lowest index is returned.

floats/floats_test.go

@@ -760,6 +760,68 @@ func TestMulTo(t *testing.T) {
 	}
 }
 
+func TestNaN(t *testing.T) {
+	tests := []struct {
+		payload uint64
+		bits    uint64
+	}{
+		{0, math.Float64bits(0 / func() float64 { return 0 }())}, // Hide the division by zero from the compiler.
+		{1, math.Float64bits(math.NaN())},
+		{1954, 0x7ff80000000007a2}, // R NA.
+	}
+
+	for _, test := range tests {
+		nan := NaN(test.payload)
+		if !math.IsNaN(nan) {
+			t.Errorf("expected NaN value, got:%f", nan)
+		}
+
+		bits := math.Float64bits(nan)
+
+		// Strip sign bit.
+		const sign = 1 << 63
+		bits &^= sign
+		test.bits &^= sign
+
+		if bits != test.bits {
+			t.Errorf("expected NaN bit representation: got:%x want:%x", bits, test.bits)
+		}
+	}
+}
+
+func TestNaNPayload(t *testing.T) {
+	tests := []struct {
+		f       float64
+		payload uint64
+		ok      bool
+	}{
+		{0 / func() float64 { return 0 }(), 0, true}, // Hide the division by zero from the compiler.
+
+		// The following two line are written explcitly to impact potential changes to math.Copysign.
+		{math.Float64frombits(math.Float64bits(math.NaN()) | (1 << 63)), 1, true},  // math.Copysign(math.NaN(), -1)
+		{math.Float64frombits(math.Float64bits(math.NaN()) &^ (1 << 63)), 1, true}, // math.Copysign(math.NaN(), 1)
+
+		{NaN(1954), 1954, true}, // R NA.
+
+		{math.Copysign(0, -1), 0, false},
+		{0, 0, false},
+		{math.Inf(-1), 0, false},
+		{math.Inf(1), 0, false},
+
+		{math.Float64frombits(0x7ff0000000000001), 0, false}, // Signalling NaN.
+	}
+
+	for _, test := range tests {
+		payload, ok := NaNPayload(test.f)
+		if payload != test.payload {
+			t.Errorf("expected NaN payload: got:%x want:%x", payload, test.payload)
+		}
+		if ok != test.ok {
+			t.Errorf("expected NaN status: got:%t want:%t", ok, test.ok)
+		}
+	}
+}
+
 func TestNearest(t *testing.T) {
 	s := []float64{6.2, 3, 5, 6.2, 8}
 	ind := Nearest(s, 2.0)