mirror of git://git.musl-libc.org/musl
math: fix scalbn and scalbnf on overflow/underflow
old code was correct only if the result was stored (without the excess precision) or musl was compiled with -ffloat-store. (see note 160 in n1570.pdf section 6.8.6.4)
This commit is contained in:
parent
662ed20065
commit
666271c105
|
@ -2,6 +2,8 @@
|
|||
|
||||
double scalbn(double x, int n)
|
||||
{
|
||||
/* make sure result is stored as double on overflow or underflow */
|
||||
volatile double z;
|
||||
double scale;
|
||||
|
||||
if (n > 1023) {
|
||||
|
@ -10,8 +12,10 @@ double scalbn(double x, int n)
|
|||
if (n > 1023) {
|
||||
x *= 0x1p1023;
|
||||
n -= 1023;
|
||||
if (n > 1023)
|
||||
return x * 0x1p1023;
|
||||
if (n > 1023) {
|
||||
z = x * 0x1p1023;
|
||||
return z;
|
||||
}
|
||||
}
|
||||
} else if (n < -1022) {
|
||||
x *= 0x1p-1022;
|
||||
|
@ -19,10 +23,13 @@ double scalbn(double x, int n)
|
|||
if (n < -1022) {
|
||||
x *= 0x1p-1022;
|
||||
n += 1022;
|
||||
if (n < -1022)
|
||||
return x * 0x1p-1022;
|
||||
if (n < -1022) {
|
||||
z = x * 0x1p-1022;
|
||||
return z;
|
||||
}
|
||||
}
|
||||
}
|
||||
INSERT_WORDS(scale, (uint32_t)(0x3ff+n)<<20, 0);
|
||||
return x * scale;
|
||||
z = x * scale;
|
||||
return z;
|
||||
}
|
||||
|
|
|
@ -2,6 +2,8 @@
|
|||
|
||||
float scalbnf(float x, int n)
|
||||
{
|
||||
/* make sure result is stored as double on overflow or underflow */
|
||||
volatile float z;
|
||||
float scale;
|
||||
|
||||
if (n > 127) {
|
||||
|
@ -10,8 +12,10 @@ float scalbnf(float x, int n)
|
|||
if (n > 127) {
|
||||
x *= 0x1p127f;
|
||||
n -= 127;
|
||||
if (n > 127)
|
||||
return x * 0x1p127f;
|
||||
if (n > 127) {
|
||||
z = x * 0x1p127f;
|
||||
return z;
|
||||
}
|
||||
}
|
||||
} else if (n < -126) {
|
||||
x *= 0x1p-126f;
|
||||
|
@ -19,10 +23,13 @@ float scalbnf(float x, int n)
|
|||
if (n < -126) {
|
||||
x *= 0x1p-126f;
|
||||
n += 126;
|
||||
if (n < -126)
|
||||
return x * 0x1p-126f;
|
||||
if (n < -126) {
|
||||
z = x * 0x1p-126f;
|
||||
return z;
|
||||
}
|
||||
}
|
||||
}
|
||||
SET_FLOAT_WORD(scale, (uint32_t)(0x7f+n)<<23);
|
||||
return x * scale;
|
||||
z = x * scale;
|
||||
return z;
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue