git: dd7dbb8021ef - stable/13 - Update scalbn* functions to the musl versions
Alex Richardson
arichardson at FreeBSD.org
Thu Apr 22 11:08:19 UTC 2021
The branch stable/13 has been updated by arichardson:
URL: https://cgit.FreeBSD.org/src/commit/?id=dd7dbb8021ef958d8a9f7417812787e2f5931da2
commit dd7dbb8021ef958d8a9f7417812787e2f5931da2
Author: Alex Richardson <arichardson at FreeBSD.org>
AuthorDate: 2021-02-24 16:56:53 +0000
Commit: Alex Richardson <arichardson at FreeBSD.org>
CommitDate: 2021-04-22 09:42:12 +0000
Update scalbn* functions to the musl versions
The only diff compared to musl is a minor change to scalbnl() to replace
musl's union ldshape with union IEEEl2bits.
This fixes the scalbn tests on non-x86 (since x86 has an assembly version
that is used instead).
Musl commit messages:
commit 8c44a060243f04283ca68dad199aab90336141db
Author: Szabolcs Nagy <nsz at port70.net>
Date: Mon Apr 3 02:38:13 2017 +0200
fix scalbn when result is in the subnormal range
in nearest rounding mode scalbn could introduce double rounding error
when an intermediate value and the final result were both in the
subnormal range e.g.
scalbn(0x1.7ffffffffffffp-1, -1073)
returned 0x1p-1073 instead of 0x1p-1074, because the intermediate
computation got rounded to 0x1.8p-1023.
with the fix an intermediate value can only be in the subnormal range
if the final result is 0 which is correct even after double rounding.
(there still can be two roundings so signals may be raised twice, but
that's only observable with trapping exceptions which is not supported.)
commit 2eaed464e2080d8321d3903b71086a1ecfc4ee4a
Author: Szabolcs Nagy <nsz at port70.net>
Date: Wed Sep 4 15:52:54 2013 +0000
math: use float_t and double_t in scalbnf and scalbn
remove STRICT_ASSIGN (c99 semantics is assumed) and use the conventional
union to prepare the scaling factor (so libm.h is no longer needed)
commit 1b77b9072f374bd26eb0574b83a0d5f18d75ec60
Author: Szabolcs Nagy <nsz at port70.net>
Date: Thu Aug 15 10:07:46 2013 +0000
math: minor scalbn*.c simplification
commit c4359e01303da2755fe7e8033826b132eb3659b1
Author: Szabolcs Nagy <nsz at port70.net>
Date: Tue Nov 13 10:55:35 2012 +0100
math: excess precision fix modf, modff, scalbn, scalbnf
old code was correct only if the result was stored (without the
excess precision) or musl was compiled with -ffloat-store.
now we use STRICT_ASSIGN to work around the issue.
(see note 160 in c11 section 6.8.6.4)
commit 666271c105e4137bdfa195e217799d74143370d4
Author: Szabolcs Nagy <nsz at port70.net>
Date: Tue Nov 13 10:30:40 2012 +0100
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)
commit 8051e08e10d2b739fcfcbc6bc7466e8d77fa49f1
Author: nsz <nsz at port70.net>
Date: Mon Mar 19 10:54:07 2012 +0100
simplify scalbn*.c implementations
The old scalbn.c was wrong and slow, the new one is just slow.
(scalbn(0x1p+1023,-2097) should give 0x1p-1074, but the old code gave 0)
Reviewed By: dim
Differential Revision: https://reviews.freebsd.org/D28872
(cherry picked from commit 760b2ffc552985017c5bc74fb9966d4dc1b5b691)
---
lib/msun/src/s_scalbn.c | 86 +++++++++++++++-----------------------------
lib/msun/src/s_scalbnf.c | 78 ++++++++++++++--------------------------
lib/msun/src/s_scalbnl.c | 92 ++++++++++++++++++------------------------------
3 files changed, 90 insertions(+), 166 deletions(-)
diff --git a/lib/msun/src/s_scalbn.c b/lib/msun/src/s_scalbn.c
index b048b0596b67..219cd8f0c989 100644
--- a/lib/msun/src/s_scalbn.c
+++ b/lib/msun/src/s_scalbn.c
@@ -1,63 +1,35 @@
-/* @(#)s_scalbn.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+#include <math.h>
+#include <stdint.h>
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-/*
- * scalbn (double x, int n)
- * scalbn(x,n) returns x* 2**n computed by exponent
- * manipulation rather than by actually performing an
- * exponentiation or a multiplication.
- */
-
-#include <float.h>
-
-#include "math.h"
-#include "math_private.h"
-
-static const double
-two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
-twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
-huge = 1.0e+300,
-tiny = 1.0e-300;
-
-double
-scalbn (double x, int n)
+double scalbn(double x, int n)
{
- int32_t k,hx,lx;
- EXTRACT_WORDS(hx,lx,x);
- k = (hx&0x7ff00000)>>20; /* extract exponent */
- if (k==0) { /* 0 or subnormal x */
- if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
- x *= two54;
- GET_HIGH_WORD(hx,x);
- k = ((hx&0x7ff00000)>>20) - 54;
- if (n< -50000) return tiny*x; /*underflow*/
- }
- if (k==0x7ff) return x+x; /* NaN or Inf */
- k = k+n;
- if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
- if (k > 0) /* normal result */
- {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
- if (k <= -54) {
- if (n > 50000) /* in case integer overflow in n+k */
- return huge*copysign(huge,x); /*overflow*/
- else
- return tiny*copysign(tiny,x); /*underflow*/
+ union {double f; uint64_t i;} u;
+ double_t y = x;
+
+ if (n > 1023) {
+ y *= 0x1p1023;
+ n -= 1023;
+ if (n > 1023) {
+ y *= 0x1p1023;
+ n -= 1023;
+ if (n > 1023)
+ n = 1023;
+ }
+ } else if (n < -1022) {
+ /* make sure final n < -53 to avoid double
+ rounding in the subnormal range */
+ y *= 0x1p-1022 * 0x1p53;
+ n += 1022 - 53;
+ if (n < -1022) {
+ y *= 0x1p-1022 * 0x1p53;
+ n += 1022 - 53;
+ if (n < -1022)
+ n = -1022;
+ }
}
- k += 54; /* subnormal result */
- SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
- return x*twom54;
+ u.i = (uint64_t)(0x3ff+n)<<52;
+ x = y * u.f;
+ return x;
}
#if (LDBL_MANT_DIG == 53)
diff --git a/lib/msun/src/s_scalbnf.c b/lib/msun/src/s_scalbnf.c
index 21d001c1faf1..3a46470b5661 100644
--- a/lib/msun/src/s_scalbnf.c
+++ b/lib/msun/src/s_scalbnf.c
@@ -1,57 +1,33 @@
-/* s_scalbnf.c -- float version of s_scalbn.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian at cygnus.com.
- */
+#include <math.h>
+#include <stdint.h>
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-#include "math.h"
-#include "math_private.h"
-
-static const float
-two25 = 3.355443200e+07, /* 0x4c000000 */
-twom25 = 2.9802322388e-08, /* 0x33000000 */
-huge = 1.0e+30,
-tiny = 1.0e-30;
-
-float
-scalbnf (float x, int n)
+float scalbnf(float x, int n)
{
- int32_t k,ix;
- GET_FLOAT_WORD(ix,x);
- k = (ix&0x7f800000)>>23; /* extract exponent */
- if (k==0) { /* 0 or subnormal x */
- if ((ix&0x7fffffff)==0) return x; /* +-0 */
- x *= two25;
- GET_FLOAT_WORD(ix,x);
- k = ((ix&0x7f800000)>>23) - 25;
- if (n< -50000) return tiny*x; /*underflow*/
- }
- if (k==0xff) return x+x; /* NaN or Inf */
- k = k+n;
- if (k > 0xfe) return huge*copysignf(huge,x); /* overflow */
- if (k > 0) /* normal result */
- {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
- if (k <= -25) {
- if (n > 50000) /* in case integer overflow in n+k */
- return huge*copysignf(huge,x); /*overflow*/
- else
- return tiny*copysignf(tiny,x); /*underflow*/
+ union {float f; uint32_t i;} u;
+ float_t y = x;
+
+ if (n > 127) {
+ y *= 0x1p127f;
+ n -= 127;
+ if (n > 127) {
+ y *= 0x1p127f;
+ n -= 127;
+ if (n > 127)
+ n = 127;
+ }
+ } else if (n < -126) {
+ y *= 0x1p-126f * 0x1p24f;
+ n += 126 - 24;
+ if (n < -126) {
+ y *= 0x1p-126f * 0x1p24f;
+ n += 126 - 24;
+ if (n < -126)
+ n = -126;
+ }
}
- k += 25; /* subnormal result */
- SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
- return x*twom25;
+ u.i = (uint32_t)(0x7f+n)<<23;
+ x = y * u.f;
+ return x;
}
__strong_reference(scalbnf, ldexpf);
diff --git a/lib/msun/src/s_scalbnl.c b/lib/msun/src/s_scalbnl.c
index 28b0cf9b7294..65a9415b0d21 100644
--- a/lib/msun/src/s_scalbnl.c
+++ b/lib/msun/src/s_scalbnl.c
@@ -1,18 +1,7 @@
-/* @(#)s_scalbn.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
+#include <math.h>
+#include <float.h>
+#include "math_private.h"
+#include "fpmath.h"
/*
* scalbnl (long double x, int n)
* scalbnl(x,n) returns x* 2**n computed by exponent
@@ -20,52 +9,39 @@ __FBSDID("$FreeBSD$");
* exponentiation or a multiplication.
*/
-/*
- * We assume that a long double has a 15-bit exponent. On systems
- * where long double is the same as double, scalbnl() is an alias
- * for scalbn(), so we don't use this routine.
- */
-
-#include <float.h>
-#include <math.h>
-
-#include "fpmath.h"
-
-#if LDBL_MAX_EXP != 0x4000
-#error "Unsupported long double format"
-#endif
-
-static const long double
-huge = 0x1p16000L,
-tiny = 0x1p-16000L;
-
-long double
-scalbnl (long double x, int n)
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double scalbnl(long double x, int n)
+{
+ return scalbn(x, n);
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double scalbnl(long double x, int n)
{
union IEEEl2bits u;
- int k;
- u.e = x;
- k = u.bits.exp; /* extract exponent */
- if (k==0) { /* 0 or subnormal x */
- if ((u.bits.manh|u.bits.manl)==0) return x; /* +-0 */
- u.e *= 0x1p+128;
- k = u.bits.exp - 128;
- if (n< -50000) return tiny*x; /*underflow*/
- }
- if (k==0x7fff) return x+x; /* NaN or Inf */
- k = k+n;
- if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow */
- if (k > 0) /* normal result */
- {u.bits.exp = k; return u.e;}
- if (k <= -128) {
- if (n > 50000) /* in case integer overflow in n+k */
- return huge*copysign(huge,x); /*overflow*/
- else
- return tiny*copysign(tiny,x); /*underflow*/
+
+ if (n > 16383) {
+ x *= 0x1p16383L;
+ n -= 16383;
+ if (n > 16383) {
+ x *= 0x1p16383L;
+ n -= 16383;
+ if (n > 16383)
+ n = 16383;
+ }
+ } else if (n < -16382) {
+ x *= 0x1p-16382L * 0x1p113L;
+ n += 16382 - 113;
+ if (n < -16382) {
+ x *= 0x1p-16382L * 0x1p113L;
+ n += 16382 - 113;
+ if (n < -16382)
+ n = -16382;
+ }
}
- k += 128; /* subnormal result */
- u.bits.exp = k;
- return u.e*0x1p-128;
+ u.e = 1.0;
+ u.xbits.expsign = 0x3fff + n;
+ return x * u.e;
}
+#endif
__strong_reference(scalbnl, ldexpl);
More information about the dev-commits-src-all
mailing list