The branch main has been updated by kib:
URL:
https://cgit.FreeBSD.org/src/commit/?id=d180086e6eae2e152e803ed6cf13775a7c006dc7
commit d180086e6eae2e152e803ed6cf13775a7c006dc7
Author: Steve Kargl <ka...@freebsd.org>
AuthorDate: 2025-08-12 04:26:29 +0000
Commit: Konstantin Belousov <k...@freebsd.org>
CommitDate: 2025-08-12 05:31:03 +0000
[libm] Fix undefined behavior of a left shifted of a signed integer
The patch fixes a few instances of left shifts on
signed integer entities. A 'static inline' helper function
'subnormal_ilogb()' has been added to math_private.h. This
function is then used e_fmod.c, s_ilogb(), and s_remquo.c.
The change in s_remquo.c has only been compile tested.
The change to e_fmod.c has been test on over 3 billion pairs
of subnormal numbers where testing included x > y and x < y
pairs. The test compared the output from fmod() with the
output from mpfr_fmod() from MPFR. There were no difference.
The change to s_ilogb() has had limited testing where its
output was compared against frexp(). In this testing, no
differences in output were detected.
PR: 288778
MFC after: 1 week
---
lib/msun/src/e_fmod.c | 30 ++++++++++++------------------
lib/msun/src/math_private.h | 21 +++++++++++++++++++++
lib/msun/src/s_ilogb.c | 13 +++++--------
lib/msun/src/s_remquo.c | 30 ++++++++++++------------------
4 files changed, 50 insertions(+), 44 deletions(-)
diff --git a/lib/msun/src/e_fmod.c b/lib/msun/src/e_fmod.c
index 77afd116c658..ced9cce33aa0 100644
--- a/lib/msun/src/e_fmod.c
+++ b/lib/msun/src/e_fmod.c
@@ -26,14 +26,14 @@ static const double one = 1.0, Zero[] = {0.0, -0.0,};
double
fmod(double x, double y)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
- u_int32_t lx,ly,lz;
+ int32_t hx, hy, hz, ix, iy, n, sx;
+ u_int32_t lx, ly, lz;
EXTRACT_WORDS(hx,lx,x);
EXTRACT_WORDS(hy,ly,y);
sx = hx&0x80000000; /* sign of x */
- hx ^=sx; /* |x| */
- hy &= 0x7fffffff; /* |y| */
+ hx ^= sx; /* |x| */
+ hy &= 0x7fffffff; /* |y| */
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
@@ -46,22 +46,16 @@ fmod(double x, double y)
}
/* determine ix = ilogb(x) */
- if(hx<0x00100000) { /* subnormal x */
- if(hx==0) {
- for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
- } else {
- for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
- }
- } else ix = (hx>>20)-1023;
+ if(hx<0x00100000)
+ ix = subnormal_ilogb(hx, lx);
+ else
+ ix = (hx>>20)-1023;
/* determine iy = ilogb(y) */
- if(hy<0x00100000) { /* subnormal y */
- if(hy==0) {
- for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
- } else {
- for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
- }
- } else iy = (hy>>20)-1023;
+ if(hy<0x00100000)
+ iy = subnormal_ilogb(hy, ly);
+ else
+ iy = (hy>>20)-1023;
/* set up {hx,lx}, {hy,ly} and align y to x */
if(ix >= -1022)
diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
index 1595f902846c..0711f2f33c41 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -739,6 +739,27 @@ irintl(long double x)
(ar) = (x) - (ai); \
} while (0)
+/*
+ * For a double entity split into high and low parts, compute ilogb.
+ */
+static inline int32_t
+subnormal_ilogb(int32_t hi, int32_t lo)
+{
+ int32_t j;
+ uint32_t i;
+
+ j = -1022;
+ if (hi == 0) {
+ j -= 21;
+ i = (uint32_t)lo;
+ } else
+ i = (uint32_t)hi << 11;
+
+ for (; i < 0x7fffffff; i <<= 1) j -= 1;
+
+ return (j);
+}
+
#ifdef DEBUG
#if defined(__amd64__) || defined(__i386__)
#define breakpoint() asm("int $3")
diff --git a/lib/msun/src/s_ilogb.c b/lib/msun/src/s_ilogb.c
index 27e0bbb8735b..aa707d51d7b7 100644
--- a/lib/msun/src/s_ilogb.c
+++ b/lib/msun/src/s_ilogb.c
@@ -21,21 +21,18 @@
#include "math.h"
#include "math_private.h"
- int ilogb(double x)
+int
+ilogb(double x)
{
- int32_t hx,lx,ix;
+ int32_t hx, ix, lx;
EXTRACT_WORDS(hx,lx,x);
hx &= 0x7fffffff;
if(hx<0x00100000) {
if((hx|lx)==0)
return FP_ILOGB0;
- else /* subnormal x */
- if(hx==0) {
- for (ix = -1043; lx>0; lx<<=1) ix -=1;
- } else {
- for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
- }
+ else
+ ix = subnormal_ilogb(hx, lx);
return ix;
}
else if (hx<0x7ff00000) return (hx>>20)-1023;
diff --git a/lib/msun/src/s_remquo.c b/lib/msun/src/s_remquo.c
index 206d2903cd86..b26b5619f3ad 100644
--- a/lib/msun/src/s_remquo.c
+++ b/lib/msun/src/s_remquo.c
@@ -4,7 +4,7 @@
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@@ -27,7 +27,7 @@ static const double Zero[] = {0.0, -0.0,};
double
remquo(double x, double y, int *quo)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
+ int32_t hx,hy,hz,ix,iy,n,sx;
u_int32_t lx,ly,lz,q,sxy;
EXTRACT_WORDS(hx,lx,x);
@@ -53,25 +53,19 @@ remquo(double x, double y, int *quo)
}
/* determine ix = ilogb(x) */
- if(hx<0x00100000) { /* subnormal x */
- if(hx==0) {
- for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
- } else {
- for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
- }
- } else ix = (hx>>20)-1023;
+ if(hx<0x00100000)
+ ix = subnormal_ilogb(hx, lx);
+ else
+ ix = (hx>>20)-1023;
/* determine iy = ilogb(y) */
- if(hy<0x00100000) { /* subnormal y */
- if(hy==0) {
- for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
- } else {
- for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
- }
- } else iy = (hy>>20)-1023;
+ if(hy<0x00100000)
+ iy = subnormal_ilogb(hy, ly);
+ else
+ iy = (hy>>20)-1023;
/* set up {hx,lx}, {hy,ly} and align y to x */
- if(ix >= -1022)
+ if(ix >= -1022)
hx = 0x00100000|(0x000fffff&hx);
else { /* subnormal x, shift x to normal */
n = -1022-ix;
@@ -83,7 +77,7 @@ remquo(double x, double y, int *quo)
lx = 0;
}
}
- if(iy >= -1022)
+ if(iy >= -1022)
hy = 0x00100000|(0x000fffff&hy);
else { /* subnormal y, shift y to normal */
n = -1022-iy;
