Ok, ok, I've fixed gradual underflow on packing too.
Still don't know what to do with Inf/NaN, so I raise an exception on pack(), and something random decided by ldexp() will happen on unpack().
This commit is contained in:
parent
6083f0e9ce
commit
4ccc531f34
|
@ -207,7 +207,7 @@ get_ulong(v, p)
|
|||
Point Arithmetic). See the following URL:
|
||||
http://www.psc.edu/general/software/packages/ieee/ieee.html */
|
||||
|
||||
/* XXX Signed zero, infinity, underflow, NaN are not handled quite right? */
|
||||
/* XXX Inf/NaN are not handled quite right (but underflow is!) */
|
||||
|
||||
static int
|
||||
pack_float(x, p, incr)
|
||||
|
@ -217,8 +217,8 @@ pack_float(x, p, incr)
|
|||
{
|
||||
int s;
|
||||
int e;
|
||||
double fl;
|
||||
long f;
|
||||
double f;
|
||||
long fbits;
|
||||
|
||||
if (x < 0) {
|
||||
s = 1;
|
||||
|
@ -226,13 +226,15 @@ pack_float(x, p, incr)
|
|||
}
|
||||
else
|
||||
s = 0;
|
||||
fl = frexp(x, &e);
|
||||
/* Normalize fl to be in the range [1.0, 2.0) */
|
||||
if (0.5 <= fl && fl < 1.0) {
|
||||
fl *= 2.0;
|
||||
|
||||
f = frexp(x, &e);
|
||||
|
||||
/* Normalize f to be in the range [1.0, 2.0) */
|
||||
if (0.5 <= f && f < 1.0) {
|
||||
f *= 2.0;
|
||||
e--;
|
||||
}
|
||||
else if (fl == 0.0) {
|
||||
else if (f == 0.0) {
|
||||
e = 0;
|
||||
}
|
||||
else {
|
||||
|
@ -240,41 +242,40 @@ pack_float(x, p, incr)
|
|||
"frexp() result out of range");
|
||||
return -1;
|
||||
}
|
||||
e += 127;
|
||||
if (e >= 255) {
|
||||
/* XXX 255 itself is reserved for Inf/NaN */
|
||||
|
||||
if (e >= 128) {
|
||||
/* XXX 128 itself is reserved for Inf/NaN */
|
||||
PyErr_SetString(PyExc_OverflowError,
|
||||
"float too large to pack with f format");
|
||||
return -1;
|
||||
}
|
||||
else if (e <= 0) {
|
||||
/* XXX Underflow -- could do better, but who cares? */
|
||||
fl = 0.0;
|
||||
else if (e < -126) {
|
||||
/* Gradual underflow */
|
||||
f = ldexp(f, 126 + e);
|
||||
e = 0;
|
||||
}
|
||||
if (fl == 0.0) {
|
||||
f = 0;
|
||||
}
|
||||
else {
|
||||
fl -= 1.0; /* Get rid of leading 1 */
|
||||
fl *= 8388608.0; /* 2**23 */
|
||||
f = (long) floor(fl + 0.5); /* Round */
|
||||
e += 127;
|
||||
f -= 1.0; /* Get rid of leading 1 */
|
||||
}
|
||||
|
||||
f *= 8388608.0; /* 2**23 */
|
||||
fbits = (long) floor(f + 0.5); /* Round */
|
||||
|
||||
/* First byte */
|
||||
*p = (s<<7) | (e>>1);
|
||||
p += incr;
|
||||
|
||||
/* Second byte */
|
||||
*p = ((e&1)<<7) | (f>>16);
|
||||
*p = ((e&1)<<7) | (fbits>>16);
|
||||
p += incr;
|
||||
|
||||
/* Third byte */
|
||||
*p = (f>>8) & 0xFF;
|
||||
*p = (fbits>>8) & 0xFF;
|
||||
p += incr;
|
||||
|
||||
/* Fourth byte */
|
||||
*p = f&0xFF;
|
||||
*p = fbits&0xFF;
|
||||
|
||||
/* Done */
|
||||
return 0;
|
||||
|
@ -288,7 +289,7 @@ pack_double(x, p, incr)
|
|||
{
|
||||
int s;
|
||||
int e;
|
||||
double fl;
|
||||
double f;
|
||||
long fhi, flo;
|
||||
|
||||
if (x < 0) {
|
||||
|
@ -297,13 +298,15 @@ pack_double(x, p, incr)
|
|||
}
|
||||
else
|
||||
s = 0;
|
||||
fl = frexp(x, &e);
|
||||
/* Normalize fl to be in the range [1.0, 2.0) */
|
||||
if (0.5 <= fl && fl < 1.0) {
|
||||
fl *= 2.0;
|
||||
|
||||
f = frexp(x, &e);
|
||||
|
||||
/* Normalize f to be in the range [1.0, 2.0) */
|
||||
if (0.5 <= f && f < 1.0) {
|
||||
f *= 2.0;
|
||||
e--;
|
||||
}
|
||||
else if (fl == 0.0) {
|
||||
else if (f == 0.0) {
|
||||
e = 0;
|
||||
}
|
||||
else {
|
||||
|
@ -311,31 +314,30 @@ pack_double(x, p, incr)
|
|||
"frexp() result out of range");
|
||||
return -1;
|
||||
}
|
||||
e += 1023;
|
||||
if (e >= 2047) {
|
||||
/* XXX 2047 itself is reserved for Inf/NaN */
|
||||
|
||||
if (e >= 1024) {
|
||||
/* XXX 1024 itself is reserved for Inf/NaN */
|
||||
PyErr_SetString(PyExc_OverflowError,
|
||||
"float too large to pack with d format");
|
||||
return -1;
|
||||
}
|
||||
else if (e <= 0) {
|
||||
/* XXX Underflow -- could do better, but who cares? */
|
||||
fl = 0.0;
|
||||
else if (e < -1022) {
|
||||
/* Gradual underflow */
|
||||
f = ldexp(f, 1022 + e);
|
||||
e = 0;
|
||||
}
|
||||
if (fl == 0.0) {
|
||||
fhi = flo = 0;
|
||||
}
|
||||
else {
|
||||
/* fhi receives the high 28 bits; flo the low 24 bits */
|
||||
fl -= 1.0;
|
||||
fl *= 268435456.0; /* 2**28 */
|
||||
fhi = (long) floor(fl); /* Truncate */
|
||||
fl -= (double)fhi;
|
||||
fl *= 16777216.0; /* 2**24 */
|
||||
flo = (long) floor(fl + 0.5); /* Round */
|
||||
e += 1023;
|
||||
f -= 1.0; /* Get rid of leading 1 */
|
||||
}
|
||||
|
||||
/* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */
|
||||
f *= 268435456.0; /* 2**28 */
|
||||
fhi = (long) floor(f); /* Truncate */
|
||||
f -= (double)fhi;
|
||||
f *= 16777216.0; /* 2**24 */
|
||||
flo = (long) floor(f + 0.5); /* Round */
|
||||
|
||||
/* First byte */
|
||||
*p = (s<<7) | (e>>4);
|
||||
p += incr;
|
||||
|
|
Loading…
Reference in New Issue