Date: Thu, 22 Apr 2021 11:08:18 GMT From: Alex Richardson <arichardson@FreeBSD.org> To: src-committers@FreeBSD.org, dev-commits-src-all@FreeBSD.org, dev-commits-src-branches@FreeBSD.org Subject: git: dd7dbb8021ef - stable/13 - Update scalbn* functions to the musl versions Message-ID: <202104221108.13MB8IFS088407@gitrepo.freebsd.org>
next in thread | raw e-mail | index | archive | help
The branch stable/13 has been updated by arichardson: URL: https://cgit.FreeBSD.org/src/commit/?id=dd7dbb8021ef958d8a9f7417812787e2f5931da2 commit dd7dbb8021ef958d8a9f7417812787e2f5931da2 Author: Alex Richardson <arichardson@FreeBSD.org> AuthorDate: 2021-02-24 16:56:53 +0000 Commit: Alex Richardson <arichardson@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@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@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@port70.net> Date: Thu Aug 15 10:07:46 2013 +0000 math: minor scalbn*.c simplification commit c4359e01303da2755fe7e8033826b132eb3659b1 Author: Szabolcs Nagy <nsz@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@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@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@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);
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?202104221108.13MB8IFS088407>