Oops! Spoke too soon.
There is a new big-endian bug:
- *Quad_RVAL.tc:* DOMAIN ERROR in 2⍟ on random 32-bit floats.
- *Domino_stress.tc:* DOMAIN ERROR in A⌹B on random matrices.
*Cause:* Quad_RVAL::random_ieee() builds a random double by writing the
bytes of a double/char[8] union, assuming little-endian layout. On sparc64
the exponent lands in the wrong bytes. Asking for 6 random floats there
gives
5.29E45 ¯1 ¯9.5E212 65329.9 ¯9.6E165 ¯2.7E143; amd64 gives values in [0,1).
*Fix:* a small patch, below. It assembles the bits in a uint64_t and
memcpy()s them into the double. It takes the same bits of the random number
as before, so the random sequence is unchanged on amd64; I checked that
with a fixed seed. On amd64 it builds and both testcases pass.
There is also a failure due to a bug in a test which hard-codes Linux
socket constants (SOL_SOCKET=1, SO_REUSEADDR=2) which fails on sparc64.
I'll patch that too when I get a chance.
--- a/trunk/src/Quad_RVAL.cc
+++ b/trunk/src/Quad_RVAL.cc
@@ -22,6 +22,7 @@
*/
#include <stdlib.h>
+#include <string.h>
#include "Avec.hh"
#include "Common.hh"
@@ -736,38 +737,34 @@
double
Quad_RVAL::random_ieee()
{
-union { double f;
- char bytes[8];
- } u;
+double f;
do {
const int64_t rand64 = rand17()
^ (rand17() << 16)
^ (rand17() << 32)
^ (rand17() << 48);
-
+
enum { BIAS = 1023, // IEEE
EXPO = 0 + BIAS, // 2^0
- EXPO_LSB = EXPO & 0x0F, // 0..15
- EXPO_MSB = EXPO >> 4 // 0..15
};
- u.bytes[7] = EXPO_MSB; // SIGN + and exponent MSBs
- u.bytes[6] = (EXPO_LSB << 4) // exponent LSBs
- | (rand64 >> 56 & 0x0F); // 52 bit fraction MSBs
- u.bytes[5] = rand64 >> 40; // 52 bit fraction
- u.bytes[4] = rand64 >> 32; // 52 bit fraction
- u.bytes[3] = rand64 >> 24; // 52 bit fraction
- u.bytes[2] = rand64 >> 16; // 52 bit fraction
- u.bytes[1] = rand64 >> 8; // 52 bit fraction
- u.bytes[0] = rand64 >> 0; // 52 bit fraction LSBs
- } while (!isnormal(u.f));
+ // sign +, exponent EXPO, and a random 52 bit fraction (bits 0-47 and
+ // 56-59 of rand64). Assemble the bits in a uint64_t and memcpy() them
+ // into f, rather than storing bytes into a union, so that the result
+ // does not depend on the host byte order.
+ //
+ const uint64_t bits = uint64_t(EXPO) << 52
+ | (uint64_t(rand64) >> 56 & 0x0F) << 48
+ | (uint64_t(rand64) & 0xFFFFFFFFFFFFULL);
+ memcpy(&f, &bits, sizeof(f));
+ } while (!isnormal(f));
- // at this point: 1.0 < u.f < 2.0
+ // at this point: 1.0 < f < 2.0
//
if (desired_purpose == 1) // APL2 mode: an IEEE 32-bit (float) value
- return double(float(u.f - 1.0));
+ return double(float(f - 1.0));
- return u.f - 1.0;
+ return f - 1.0;
}
//────────────────────────────────────────────────────────────────────────────
void
